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ABSTRACT 


In this study, the basic Muskingum method, the Musk ingum -Cunge 
method CMC) and the various improvements of MC method are 
discussed. At first, MC method of flood routing is compared with 
that of the implicit Preissmann scheme. It has been found that the 
MC method has less attenuation and more translation of outflow peak 
discharge. The MC method predicts more initial flows than the 
implicit Preissmann scheme. 

To improve the MC method, an attempt is made by changing its 
forward finite difference scheme into a Leap frog finite difference 


scheme, named 

MCLF model. 

It 

has been 

found 

that 

this 

proposed 

model is 

of 

second order 

accurate 

scheme 

and 

the 

numer ical 

diffusion 

is 

eliminated 

considerably 

. Further, 

the 

solution 


approaches the analytical solution of the kinematic wave equation, 
creating very less attenuation of peak discharge and hence do not 
adequately represent the physical phenomena. 

To predict the better initial flows by the MC method, a method 
of segregating IH which uses the MC method called SMC modal is 
proposed- The SMC model predicts the flows upto the first Q more 

- !50 

accurately than the MC method. 



iv 


A different apprciach is otade to improve the MC method by using 
the Laplace transform solution of a step input function for 
Muskingum equation with the relevant boundary conditions. This 
model is named as MLT model. In the MLT model, the given IH is 
approximated to a input staircase function and routed explicitly. 
The MLT method predicts outflow peaks to the same degree of 
accuracy as the MC method. The time of occurrence of the routed 

peak outflow is also predicted to the same degree of accuracy as in 

> 

the MC method. The MLT method is thus a viable alternative to the 
MC method. 
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NOTATION 


cross-sectional area 
breadth of channel 
celerity 
Courant number 
acceleration due to gravity 
distance node 

Inflow ordinate at time t = 0 
inflow at any time 
time node 

average reach travel time 
length of reach 

Manning's roughness coefficient 

variable to transform time domain to Laplace domain 

outflow at any time 

peak discharge 

Storage 

bed slope 

t ime 

flow velocity 

distance in direction of flow 
depth of flow 
time increment 
distance increment 


Coefficient reflecting effect of Inflow and outflow 



CHAPTER-1 


INTRODUCTION 

Flood routing is the technique of determining the flood 
hydrograph at a section of river by utilizing the data of flood 
flow at one or more upstream sections. The hydrologic analysis of 
problems such as flood forecasting, flood protection, reservoir 
design and spillway design invariably include flood routing. In 
these applications, two broad categories of routing can be 
recognized. These are: 

1. Reservoir routing. 

2. Channel routing. 

In reservoir routing the effect of flood wave entering a 

reservoir is studied. In channel routing the changes in the shape 

/ 

of a hydrograph as it travels down a channel is studied. Flood 
routing for channels may be done by using the systems approach, 
also called as hydrologic flood routing. If flood routing is done 
by a process approach, it is called hydraulic flood routing. The 
hydrologic modeling involves a conceptual modeling whereas the 
later is done by mathematical modeling. Hydrologic channel routing 
is based on the storage concept. The hydraulic models are based on 
the solution of a pair of non-linear partial differential 
equations, the Saint Venant equations, based on the principles of 
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also called as hydrologic flood routing. If flood routing is done 

I 

by a process approach, it is called hydraulic flood routing. The 
hydrologic modeling involves a conceptual modeling whereas the 
later is done by mathematical modeling. Hydrologic channel routing 
is based on the storage concept. The hydraulic otodels are based on 
the solution of a pair of non-linear partial differential 
equations, the Saint Venant equations, based on the principles of 



mass and momentum conservation. 


2 


There are broadly two techniques to solve the hydrodynamic 

Saint -Venant equations 

1. Approximate methods which are based on equations of continuity 
and on a drastically curtailed equations of motion- 

2. Complete Numerical methods which solve the basic Saint-Venant 
equations through numerical model ingC233 . 

Thfeic: * hybrid approach which is similar in nature to hyrologic 
■outing methods, yet contains sufficient physical information to 
lompare favorably with the more complex hydraulic routing 
echniques. This hybrid approach is the basis of lluskinguro-Cunge 
lethod of flood routing. 

The Muskingum method is a hydrologic model which is more a 
onceptual model with empirically determined storage functions. The 
irect relationship between the steady .flow rating function 
= FCyJ and the storage function provides a basis for this mode. 

The Muskingum-Cunge method is a hybrid model in which the 
luting parameters are calculated as a function of 
ysical properties: 

1) Reach, length Ax 

2) Peak diechax-ge Q 

y p 

3) KirxBmntic wax>B cBl&rity c 

Bottom slopo s . 

o 


the following 



The various flood routing methods are shown in the Flow 


hart -1. 


Flood routing 


Channel routing 


Reservoir routing 

1. Modified Puls method, 

2. GoodRich method . 


Hydrologic routing 
1. Musk ingum method 


Hybrid routing 
M-Cunge method 


Hydraulic routing 


Method of Finite difference Finite element Boundary Spectral 
Character sties methods methods integral methods 

11 methods 


r 


"1 


Explicit methods 
l.MacCormack scheme 
2- Lambda Scheme 
S.Gabutti scheme 

a,3 ate examples 

Flow Chart for methods of flood routing 


Implicit methods 

1. Preismann Scheme. 

2, Abott^'s Scheiinc. 


The basic Muskingum method, the Masking urn -Cunge method and its 
various variations t improvements) are discussed in Chapter 2. In 
Chapter analytical method which is developed by using a 
Leap-frog scheme is presented. A technique of Muskingum -Cunge 
routing by segregating the input hydrograph is developed and its 
character ist ics are discussed in CJiapter 4. An approach for routing 
the flood by a lumped model is empirically developed by making use 
of the solution to a step input function of Muskingum equation with 
suitable boundary conditions in Chapter 5. Chapter 6 discusses the 
significant conclusions of this study. 




CHAPTER-2 


LITERATURE REVIEW 


2.0 INTRODUCTION: 

A general form of the resistance formulas, such as Chezy or 
Manning may be written as 

Q = C' A R™ / s^ C2.1) 

where C' is the empirical constant, R is the hydraulic radius, m is 
an exponent and s^. is the slope of the energy grade line. In 
Chezy's formula , C' = C, m = 0.5 . For steady uniform flows, Q is 
the normal discharge Q and s, = s . 

n f O 

2.1 KINEMATIC WAVE MODEL: 


For steady uniform flows EqC2. 1) becomes 
Q = C' A r”’ VsT^ 

n O 

Eliminating C''AR*" from EqC2.1) and Eq(2.2) , we obtain 

Q y 

Q = !1_ Vs/ 

/V ' 

From the dynamic equation 

dy 


s, = s - — 3 — 
f o dx 


V dV 
g dx 


dv 


g dt 

Substitution of expression for s in Eq(2.31 yields 


C2.2) 


C2.31 


(2.4) 


Q = Q / 1 I L_ ^ 

n S dx 


V dv" 

s g dx 
o 


T Bv 


s g dt 
o 


Kinematic 

Diffusion 


J'^'nanrrilc 




(2.5) 



The Eq(2.5) demonstrates the distinguishing characteristic of 
the kinematic wave model, that the (Eg A1--2) discharge is always 
equal to the normal discharge and is thus a single valued function 


of depth. 


The St. Venant equationsiEqs. 

dA ^ ao _ 
at ^ ~ ^ 


Al-2> reduce to 

C2.6) 


and Q = Q = C A R*" -/ s^ C2.7) 

n O 

This system can be combined into a single equation .commonly 
referred to as kinematic wave equation, as 


1 

c at 


ao 

ax 


= q 


C2.8) 


in which the coefficient 1/c is either a constant (linear models) c 
a variable (non linear models). The^ parameter c, describing the 
travel speed of flood wave is called the kinematic wave speed, and 
may, at a particular cross-section and for a given Q , can be 


evaluated from Kleitz- Seddon-'s law as 
^ • dA X B dy X 


(2.9) 


2.2 MUSKINGUM-METHOD OF FLOOD ROUTING; 

In the Muskingum method of routing, the form of continuity 
equation used in the lumped model states that, for a given time 
interval, the difference between the inflow to a reach and the 
outflow from the reach must be equal to the change of storage 
within the reach. Mathematically, the storage equation is 

= I(t) - D(t) (2.10) 

, dt 

and storage - discharge relationship is given by a linear 
relat ionship. 


s = K ce I + (1- e)Q3 


( 2 . 11 ) 


where K and 0 are routing parameters. Physically, K is considered 


as the average reach travel time in the reach , and 0 is a 
coefficient reflecting the relative effects of inflow and outflow 
on each storage. 

Expressing Eqt2. 101 in the finite-difference form and 
then substituting EqC2.11} 


I + I * 0 + Q , 

, At - 


= K 




-Q^)l (2. 


12 ) 


in which At is routing period or descretization time 

interval. 

Rearranging and solving for Q .. 

+Ai 

(2.13) 


Q* A. = C I . + C I + C Q 

o i-fAt 1 i 2 i 


C = - 
o 


K 0 + 0.5 At 


(2.14a) 


C = 
1 


K 0 + 0.5 At 


(2. 14b) 


C = 
2 


K (1-0) -0.5 At 


(2. 14c) 


where 


C = K(l-0) + 0.5 At 

3 


(2.14d) 


The coefficients C , C and C are such that C + C 

0 12 o 1 


C = 1.0 
2 


2.2.0 ESTIMATION OF PARAMETERS: 

There are a number of methods available for estimating the 
parameters K and & of the Muskingum method. The available techniques 
are 

1) Graphical estimation 



2) Moment matching technique 


3) Optimization technique 
4> Least squares estimation. 

2.2.1 GRAPHICAL ESIMATION: 

In graphical estimation , a trial value of & is chosen, values 
of storage S at any time t are plotted against the corresponding 
C&I + C1-&)Q1 values. If the value of 6 is chosen correctly, a 
straight line relationship will be resulting. If an incorrect value 
of 0 is used , the plotted points will trace a looping curve. By 

trial and error, a value of 0 is so chosen that the data very 

nearly describe a straight line. The inverse of this straight line 

will give the value of K. Details are available in 

SubramanyaC243 . 


2.2.2 METHOD OF MOMENTS: 

This method is based on first determining the moments of the 
Instantaneous unit hydrograph (lUH), h(t) of the reach and then 
expressing these moments ,in terms of the moments of inflow and 
outflow . The parameters K and 0 are then estimated as €221: 


h(tl=- 


-tXtK<l-0> 


6(tl 


KCl-0)‘ 


where htt) is the lUH of the reach 


<5(t) is the Dirac <5 function or impulse inflow. 


K = M° (Q) - M® Cl) 

1 1 

{ M CO) - M Cl) 

1 

In this M^CQ), M°CI) are the first moments of 



C2. 15) 

C2. 16) 

the outflow and 



inflow hydrographs about the origin. 

M (Q) r the second moments of outflow and inflow 

2 2 

hydrographs about their respective centre of gravity. 


2.2.3 OPTIMIZATION TECHNIQUE: 


The direct optimization method C221 s based on directly 

determining the coefficients C , C , C in EqC2. 13) without first 

0 12 

estimating K and © in Eq(2. 14). This would involve minimizing the 

error function in comparing observed hydrograph with computed 

hydrograph. Again, the error function can be defined in a least 

squares sense or differently. There are in fact only two unknowns, 

since the third can be obtained from C + C +• C = 1.0. Let us 

0 12 

consider and to be the unknowns. Then Eq(2.13) reduces to 


Q* A." I* A.= C (I, .. - I 3 + C (I . - Q 3 
t+Al t+At 4 t+At t 2 t+Al t 


C2.173 


The coefficients and can be derived by comparing flows of a 
known inflow hydrograph with flows of the corresponding outflow 
hydrograph. 

Let us define 


R s; D — T 

t+At l+At 
^ = ^+Ai - ^ 

® = ^+At - °t 

Eq(2.173 can be rewritten as 

R = C F + C 6 
1 2 

which is identical in form to Eq(2. 13) with C = 0.0. 

o 


Therefore, 


ZRFSe -ERGEFG 

C = < 3 

^ Z F® S G* - C E FG 3* 


C = ( 
2 


ZReZF -SRFZFG 
o o 

Z F® Z G® - ( Z F6 3* 


Z F® Z G® 


3 



Where R is observed R. Then C 
o o 

It can be obtained 


e 

K 


C 

1 


+ 0.5 C - 
z 


TT 

i 


mr 

2 


At (C + C ) 


± 2 

T=C— 

2 


can be computed from C = 1-C -C . 

o - » z 


0.5 


C2.18) 


(2.19) 


2.2.4 LEAST SQUARES METHOD; 

The least squares method is based on minimizing the sum of 
deviations between observed storage and estimated storage for a 
given inflow and outflow sequence. That is, 

N 

E = EC S^(j) - S (j) 3^ ^ min (2.20) 

o • 

where E is the error function to be minimized, S (j) is the 

o 

observed storage for the time interval, S <j) is the 

estimated storage for the time interval, and N i* the number 
of observations or data. 

The estimated storage should be given by Eq(2.11) can be 
modified to incorporate the difference C between relative and 
absolute storage. 

S * K C0 I + (1-0) Q3+C (2.21) 

If the flood wave is starting on a dry river bed, the 
parameter C will vanish. First considering the case c 0. 
Substituting Eq(2.2i) into Eq(2.22) and dropping the index for 
brevity, 

N 

E = SCS-K0I - K(l-0) Q - C3* ^ min (2.22) 

o 

Differentiating Eq(2.22) with respect to A,B,C equating each 


time to zero, and re-arranging 



r S Q = A ZIQ + B ZQ* + CZQ 

O 

Z S = A ZI + B ZQ + CN 
o 

Solving for Ar B,C 


A 

C 


B = 




y 2 
2 1 


2 

s s 


y 2 

a 2 


y 2 

2 a 


- A ZI - B ZQ 


N 




Z S - Z I 
o 

N 



z = 
z 


2 = 
c 


Z I - 

Z IQ 

Z S G 
o 

Z IQ 
Z Q* 


t ZD* 

N 

Z Q Z I 
N 

Z S Z Q 
_ o 

N 

,Z I Z Q 

^ R — " 

CZQ)* 

N 


) 


Therefore, K= A+B 



Equetions(2.23} and .(2.24> give the parameters 
a known inflow-outflow sequence. 

If C = 0, then the error function E follows: 

N 

E = Z CS - K e I - <l-e) KQl* •> min 
1 o 

Using A and B as defined before, differentiating E 
to A and B and equating each time to zero. 

Z S I = A Z I* + B Z IQ , 
o 

ZSQ = AZIQ + BZQ* 
o 


(2.23) 


(2.24a) 

(2.24b) 

K, 9 and C for 

(2.25) 

with respect 



Solving for A and B 


11 


SIEQ*-ZSQZIQ 

A = (2.26) 

ZI* ZQ* - CZ IQ)* 

ZSQZI*-ZSIZia 

B = H 2 (2.27) 

ZI* Z Q* - (Z IQ)* 

N, 

A 

Therefore/ K=A + B / 0 = — pr — . 

COMMENT: Of all the above four methodsr the moment matching 
technique is easily adoptable and simple to use. In this study, the 
moment matching technique is used, to estimate the routing 
parameters K and © when the inflow and outflow hydrographs are 
given. 

2.3 MUSKINGUM-CUNGE METHOD: 

In the Muskingum-CungeC63 method (MC method), for the 
kinematic wave modela-^i explicit finite difference scheme had been 
proposed. The space time grid is as follows : 



Distance X 

Fig. 2.1 x-t Plans: Finite difference grid for 
Musk ing urn -Cunge Model 

In the above Fig. 2.1, the scheme introduces a weighing 
coefficient 0 , in the finite difference r epr esentat ion of the time 



derivative. Evaluating the finite difference quotients for the 


space and time derivatives at the point i+1-0, j+1/2 , the 
kinematic wave equation for the zero lateral inflowr may be written 
in finite difference form as 


, 0C I - I ) + (1-©) (Q - Q > 
^ 1 ^ Z 1 2 1 


Q + Q - I - I 

1 Z 1 z 


= 0.0 


At 2 Ax 

(2.28) 

in which <l/c> ^ f(t); and <c> is an average value of c for the 
reach. 


Introducing the travel time parameters K = 


A X 
TEJ 


the 


solution of Eq(2.28} for Q can be reduced to form identical to 

z 


the classical Muskingum coefficient equation. 


Q = CI + CI + CQ 

2 11 z z a 1 


(2. 29a) 


where , 


^1 = 




K<l-©)+ 


At 


- K © + 


At 




K(l-©)+ 


At 


(2.29b) 


(2.29c) 


K(l-©) - 


At 


C = 

3 


K(l-©) + 


At 


(2.29d) 


^ CungeCAl showed that the Muskingum scheme is a second -ord? 
approximation of the diffusion equation if the parameters are 

1 

evaluated as 

K = (2.30) 


© = 0.5 X (1- 


Q 


B s c A X 
o 


(2-31) 


and 0 <0 < 0.5. 

This Eq(2.31) for the weighting coefficient © was obtained 



from the Teylor series expansion of QCx.,yJ in the finite 

j i 

difference form of Eq(2.28) and a comparison with the coefficients 
of the diffusion equations for regular channels. 


2.3.2. Raudikivi^'s Algorithm For Estimation Of Parameters: 

Raudikivi C203 has given a simplified technique for routing of 
a flood using MC method. The algorithm for calculating the values 
of K and 9 as suggested by Raudikivi is as follows : 

Given the reach length Lr peak discharge Q whose time of 

p 

occurrence is T , breadth of channel B, whose bed slope s , 
p o 

1. Determine or assume L/T , the wave celerity c. 

p 

2. Select the sub reach length Ax. 

3. Calculate a = — i— — . 

p 2 B s 

o 

4. Calculate the radius of curvature at peak 

d*Q Q+Q-2Q 

p *■ p_ 

dt* CAt)* 

where Q ^ and arc the discharges on either side of the peak 
discharges and At is the time period of given inflow 
hydrograph. 


oi 


5. Calculate Q = 


<L/T ) 
p 


- Q 

3 p 


d^ Q 


dt' 


If Q < 0.1 Q ; then c = L/T 

p p 

If Q* > 0.1 Q redefine Q* as 

p 


Q = Q 

new p 


1 1 if H 

r-" 

o 

1 ± ' 

J 



14 


6. Q = Q - 0.5 Q* if to* < 0.1 Q ) 

P P P 

Q = Q - 0.5 Q* if (Q* > 0. 1 Q 1 

P P ^•v ^ p 

7. The parameters K and & are calculated as 

a Q 


K = 


Ax 

<c> 


and 6 = 


L Ax <c> 


Ax 

where ■jrr »nu5t be less than the value given by the curve in 

c At 

Fig. 2.2 which can be approximated as 
If 9 is less than 0.10 
y = 1422.55 0® - 283.10 0* + 20.08 0 
If 0 is greater than are equal to 0.10 

y = - 1.665 0* + 2 0 + 0.417 

. Ax 

where y = 

8. Estimate C , C and C by Eqns (2.18) 

12 3 


2.4 MUSKINBUM-CUNBE METHOD WITH VARIABLE PARAMETERS: 

Ponce and Yevjevich C183 have modified the MC method by 
defining the Courant number and cell Reynolds number in the 
following way. 

\ 

The equation for Q is as given below: 

2 

Q = CI + CI + CQ 

2 11 2 2 3 1 

Qn re-arranging the coefficients such that C = 

Q / 

which is the Courant number and D = ( b ' s ~' / ^ ^ type of cell 

o 

Reynolds number. They obtained the coefficients in terms of C and D 
as 

1 1+ C + D 


r ^ -1+ C+ D 
2 " 1+ C+ D 


(2.32b) 
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C = + D 

a 1+ C+ D 


(2.32c) 


Usually, At is fixed, and Ax and s^ are specified for each 
computational cell consisting of grid points. 

n»l n»l 



Fig. 2.3 Computational cell 

Therefore, it is necessary to determine the flood wave 
celerity c and width discharge, q for each computational cell. The 
values of c and q at grid point (j,n) are defined by 

c = 1. (2.33a) 

dA ■ j,n 

q = -i- I. (2.33b) 

B *j,r> 

in which Q = discharge; A= Flow area; B= top width. The values of C 
and D are calculated from which the routing coefficients , C^, 

C are calculated. 

3 

COMMENT: This method is a slight modification of MC method in which 
the routing parameters are estimated for each computational cell. 
This way , it becomes cumbersome, to calculate the routing 
coefficients for each cell. The data at each grid point may not be 
available. Further, it loses the simplicity of MC method. 
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2.5 SIMPLIFIED MUSKINGUM ROUTING EQUATION: 

PonceC19] modified the MC method to make the calculations 
simple. The calculations for the lumped parameter 
mode can be can be proceeded in one of two ways. The conventional 
way is to specify space and time intervals Ax and At? parameters K 
and e are calculated therefrom. An alternate way is to specify K 
and 9 and to calculate At and Ax from them. The proper choice of 
parameters can lead to a considerable 'simplification in the routing 
equation. For calculating the routing coefficients of MC if At/K 
and 6 are specified as 


At 

K~ ~ 
and © = 

It follows that C = 

* 1 

Leading to the 
equation as follows: 


1.0 
0 ; 

C = C = 1/3. 

Z 3 

simplified 


form 


of 


t2.34) 

C2.35) 

the Muskingum routing 

C2.3A) 


COMMENT; With this method it is not possible to calculate at the 
desired subreach. When one uses the a computer it loses its 
significance of easy calculations. 


2.6 MUSKINGUM METHOD AS PROPOSED BY KOUSSIS: 

KoussisE133 has partially modified the MC method. In the 
continuity equation he has carried out discretization in space only 
and retaining continuous functions for the time derivatives- Using 





the time derivative the weighted scheme indicated in Fig. 2.1, the 
above kinematic wave equation, may be written as an 
differential equation as follows: 


^ . 1 o - 1 

di Cl/cl Ax Cl-©) Cl/c] Ax Cl-©) 


I - ^ £1. (.2 37) 
^ 1-© at 


with the same assumption for <l/c> as for Eq(2.28) and assuming 
linear variation of Q over the time interval. At , the solution may 
be written in a form convenient for computation, i.e. as a 
coefficient equation of the following Muskingum type: 

Q = CI+CI + CQ 

Z It Z 2 3 1 


where the coefficients are now 

C = 


* Cc: At 


- ft ; 


(2.38a) 




(2.38b) 


C = /? 


(2.38c) 


where 


„ r cci At j 

[■ Ax(i-e) 


(2.33d> 


Koussis demonstrated that this model is a second -order 

approximation to the diffusion analogy model if 


© = 1 - 


X + 1 - r 


; and 0 < © < 0. 


(2.39) 


where 


r = 


c At 
Ax 


and X = 


B S c Ax 
o 


(2.40) 
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2.7 LEAST SQUARES PARAMETER ESTIMATION FOR MUSKINGUM FLOOD ROUTING: 

Ald<9maC13 has examined the available least square parameter 
estimation techniques for the Muskingum routing. Gill involves the 
direct use of Eq(2.11} in estimation of K and 9. Absolute values of 
storage are rarely available in practice. Gill modified Eq(2-ll) by 
setting S=S-c5';a = K©? ft - K(l-0) t where S is relative 

r r 

storage and - ct may be interpreted as the initial storage. Hence, 
the storage - discharge relation may be written as S = <? + a I .+ ft 

j 

Q^, where the circumflex denotes estimated. Accordingly, the sum of 
the squares of the relative storage estimation errors becomes 

N 2 N 2 

E = X CS S ) =.5: ( o + al.+^JQ. - S.) (2.41) 

« j=i rj rj j = l J J rj 

where S^.= the actual value of relative storage at time j At? and 


N = the total number of time steps. The minimization of E leads 


to the following explicit expressions for a and ft ; 

N N N N 

( Z I.Q. XQ. -.Z I. Z Q® ) .Z S . 

J =1 J J jsl J J = 1 J J jaj, rj 


« = D 


-1 


N N N 

+ CN.Z Q* - C .Z Q. )*3' X IS . 

J = 1 J J j=l 4 rj 

N N N N 

+ iXI. XQ.-N .ZI.Q.) ,Z Q. S.)}- 

} J=i J J J J=1 J r4 J 


(2.42) 


and 


j N N N N N 

^ = D" -^ ( .Z I. X I.Q.- .Z I* .Z Q.) .Z S . 

J J 4 4=^1 J J 4=1 rj 


N 


N 


N 


N 


+ C .Z I. .Z Q, - N .1 I.Q.) .£ I.S 

4 = a 4 4=i J 4«A J 4 4=^ 4 *^4 

+ CN .Z I* “ ( .Z I. .Z Q. S . 

J=i J J=i J J*l. J rj 


(2.43) 



where 
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N N 

N 

N N 

N 

N C .2 I® X Q* 

J = 1 J 3 

- iX I. 

J 

Q.)®1 +2 .2 I . .2 Q. 

3 J J 

.2 I.Q. 

4=1 4 4 

N M 

N 

N 


< X I .)* X Q* 

j=i J j=i i 

- .z I* 

J 

< .2 Q.)* 

4*1 J 

<2.44) 


once a and ft are computed, K and & can be determined as 

K = a + ^; © = STTT? <2.453 

but the estimated outflow, through recursive use of Eq<2. 13) as 

j —1 k 


Q. = I. 

j+i j+i 


kio 


C tl. -I. , ) 

3 j-k j-k-»-i 


where C = C + C (2.46) 

Thus, minimization of the sum of squares of the outflow rate 


estimation error 

leads to 





N j-i 

< .2 F. ,2 C 

j = l j+l k=0 2 

G ) [ .2 

4k [ j=l 

< ^2 o'* G.^) 

k = o 2 jk 

4-1 

< ,2 K 

ksO 

c’'”^ G ., ; 

2 jk 

’] 

N j-1 

- < .2 F. ^2 K 

j = l j+l k = 0 

c*'"^ e.. ) 

2 jk 

r N j-i 

.2 < , 2 C 

J»4 k*0 2 

- 

0 

■ 

CM 


and 


C 


3 


H 

X 

N 

X 


j-4 

F. ,2 C 

j+l k=0 2 


k 

( S C G 
k=0 ^2 "jk 



(2.48) 


where 

and 


F. = I.- Q. 

3 3 3 

G = I - I 

jk j-k j-k+l 


<2.49) 

<2.50) 


Eq<2.47) can be solved iteratively for C and then C can 

2 3 

computed. K and Q may be obtained in terms of C and C as 

2 3 


C At 

^ = ri-c~y 

2 


0 = 1 - 


1 + c 

a 


<2.51) 

<2.52) 


be 


COMMENT; To route a flood with the above method, considerable 
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computation effort is involved. Further, the method loses the 
simplicity of Muskingum and Muskingum-Cunge methods. 


2.8 Muskingum Method With Variable Parameters of V.P.SinghC251s 

In this model, the routing parameters of the Muskingum-^Cunge 
equation are evaluated at every discrete point on x-t plane. 
Considering on Ci+l)*”'' period, let Ki+l) be the inflow, Q(i+1) the 
outflow, wti+1) the storage, and K(i+ll the corresponding travel 
time fdr that period. 


Kti+1) = 


w(i+lJ 


CQti+l>3 


at 2m+i >ya*oo 


(2.53) 


where 




^{1-2WCX )/<4+0(> - <2m+l>/r2<l.+0(>J 

o A 

gTn pa<2Tn4.i>/<l+a> 


(2.54) 


For the Chezy formula a = 0.5 


S = bottom slope 
o 

S = water surface slope 

r- Au. - . .th period 

For the period i 


K(i) = 


w(i) 


,at 2 n)-n.>/(i+a> 


(2.55) 


[ 0 ( 1)3 

By dividing Eq(2-53) by Eq(2.5S) the following is obtained; 


K(i+1) 

_ w(i+l) 

r Qti) 

Tatznt+oxti+oo , - . . » 

I w(i+l) 

r Q(i) 1 

K(i) 

w(i) 

[ Qti+1) 

J w(i) 

[ Q(i+1) J 


(2.56) 


where 7' = Ca(2m+1) 3/(l+a) . From the kinematic wave approximation. 


Q = a A' 
z 


ftz 


Then 



Q(i+1) 

r "e A(i+1)-|^^ 

_ f L A(i+1)T^2 

r w(i+i)i 

Q(i) 

L A(i)J 

[ L A(i) J 

L J 

By substituting Eq(2.57) 

into Eq(2.56) 


K(i+1) _ 

w(i+l) r w(i) 

-,auzm+x>ft /<i+oo 


K(i) 

w(i) j_ w(i+l) 

J 
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r wCi+1) 

[ wTi 7 “ J 


C2.58J 


Dimensionless weighting coefficient 

For the Muskingum method, © can be represented as 


e - 


1 

TT“ 


(2.591 


where 1 is the characteristic length which can be epressed as 


Q 


1 = 


©H 

asr 


Q 


) = 


o 

3 

o 


d(H - H 1 
o 

— 55 


(2.60) 


where H is the bottom elevation, H is the water level, and Q 
o o 

is the flow discharge. 

For steady flow, the relationship between the water level and 
the discharge is 

H - H = h = b„ (2.61) 

o o o 

where h is the water depth, b^ is the constant, and ft is an 
exponent. Taking the derivative of Eq(2.60) 


a(H - H ) 
o 

dQ 


C b^ 3 = b ft 
oQ o o o ' o 

o 


(2.62) 


and substituting into Eq(2.59) yield 


Q 


1 = 


a(H - H ) 

O 

gg- 


G 


So ^ 


ft 




(2.63) 


Substitution of Eq(2.63) into Eq(2.59) yields 
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e = 


1 h ft gP 

-ST - “!2r-S — 


<2.64) 


The relationship between steady -flow discharge Q and 

o 

unsteady-flow discharge can be expressed as 

S, 


Q = Q < 1 + -= 
o S 


A 


<2.65) 


In general, » S. , therefore, Q Q and Eq<2.64) becomes 
O A o 


^ ft 

e = 0.5 - 


= 0.5 - Q' 


<2. 66a) 


when L = Ax, h = b^Q' , and q/c - ft h Eq<2.66a) becomes 


e = 0.5<1- 


C Ax S 


<2. 66b) 


From Eq<2.66a), the following formula can be obtained: 

^ ft 

0<i+l) = 0.5 - 2i_ S Ct3ti+1)3' 

“o ^ ° ft 

0<i) = 0.5 - CQ<i)3^ 

o 


<2. 66c) 

<2.66d) 


On substracting Eq<2.66c) from Eq<2.66d) 

0<i+l) = ©<i) + ^ ° ^ C CQ<i)3^ - CQ<i+l)3^} 

o 


COMMENT: Routing with this model is not possible unless we know 
K<1),X<1), yr OL T /9,w<l)- The field data scarcity may prevent the 

X 

use of this particular model. 

\ 

In the absence of outflow data for the catchment where inflow 
data and all the physical characteristics are known, to determine 
the accuracy of a routed flood by using any method a bench mark is 
needed. In this study, as the outflow data was not available, first 
a benchmark outflow hydrograph of the reach was created by using 
implicit Priessmann schcmeCAppendix-11. The Priessroann scheme gives 
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the numerical solution of complete Saint -Venant equations. The 
relative accuracies of the Musk inq urn -Cunqe and various other 
proposed models are studied against this kind of bench mark in 
Chapter 4 and Chapter 5. 

The implicit Priessmann scheme is first modeled and tested 
against the data sets 1 and 2 CCooley,MoinC5:3. The outflow 
hydrograph obtained from the Preissmann schemer the outflow 
hydrograph of Cooley and Moin are plotted in Figs. 2.4 and 2.5 
along with the corresponding inflow hydrograph. The numerical error 
is calculated as 

E - E 

7. error = — g — 111 XlOO 

Tn 

where E represents the peak discharge or the depth of flow at peak 
discharge. E represents the peak discharge obtained by Cooley and 

m ' 

Moin .The errors for these sets of data are shown in the Table 2.1 


Table 2.1 


Data 

set 

Typa 

Moin 

Preissmann 

% ERROR 

Set 

# 1 

Qpeak 

1764.71 

1823.24 

+ 3.33 . 

Set 

tt 2 

Peak 

Height 

29.90 

29.40 

- 1.67 

From the 

above 

Table 

2.1, as 

the error is in 

the permissible 


limits tof < 5 % 1 , the Preissmann scheme modeled is considered to 
be adequate to generate the bench mark out flow hydrographs. 



Discharge (m /s) 


! 

I 


0.0 



6«00 12.00 18.00 24.00 30.00 38.00 42.00 40.00 

Time<min) 

f’i9- 2.4 Comparison of results from Preissmann 
and Moin - Triangular IH 


30.000 


27.428 


22.88? 

0 

\ 

cn 

18.208 

e 

0 

0 ) 13.714 

d 

X 

o 

^ 8.143 




‘a.st £4.03 3e»c» 4 e.ae asi.oa ra.aa e4.a@ aa.est tm. em ise. esa i-'h.sib' 

TitneChrs) 

Fig. 2.5 Comparison of results frow Preissmann 
and Hoin - Log Pearson type III IH 



CHAPTER-3 


ANALYSIS OF KINEMATIC ROUTING 

INTRODUCTION: 

In the Muskingum-Cunge(MC) routing the partial differential 
equation for continuity is approximated by implicit 
finite-difference scheme. In this chapter the continuity equation 
is approximated by Leap Frog Finite difference scheme towards 
achieving more accuracy in explicit routing of the flood. The order 
of accuracy of the leap-frog scheme is calculated and the stability 
analysis is also investigated. The proposed model is compared with 
the MC method. 


3. 1 FIRST ORDER DIFFERENTIAL EQUATION FOR ONE-TO-ONE STAGE 
DISCHARGE RELATIONSHIP : 

The complete equations for gradually varying transient flow 
(Saint -Vcnant) equation for a rectangular prismatic channel are: 


dy 1 £Q 
dx 


= 0.0 




where 


r 

^ dv 


- Q 

1 Q| 

1 dt 

dx J 

= s - 

o &X 


X 

= distance 

along 

the river 


ytx,t) 

= depth of 

water 



(3. 1) 
(3.2) 


v(x,t) = mean velocity 



width 
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s = bed slope 

O 

QCx,t) = resistance coefficient 


1 

r + V 

av 1 

g 

1 at 

ax 1 


= slope due to inertia of water 


Eqs. (3.1) and (3.2) are two non-linear partial 
differential equations are hyperbolic and their solution is subject 
to the following conditions: 

- Two boundary conditions i.e. one upstream and one downstream of 
the considered river section. 

- An initial state at time t = t at each point of the considered 

o . ■ 

section, consisting of a depth of water y(x,t) and a velocity v(x,to) 
at each cross section. 

Assuming a single-valued stage/discharge curves for given points in 

a river, so that x = x at a given point, then: 

o 

Q| = Q(y), and hence the following relationship 

> X=Xo 

wetted area A = ACy) 

dy _ dy dA (3.3) 

A| = M-Q) (3.4) 

■x = X 

o 

dA f cJA *) £Q 

at t dQ Jxq at 

ay _ i f ^ 1 ^ 
aT “ B I dQ J at 

Substituting Eq(3.6) into Eq (3.1), we obtain the 


(3.5) 

(3.6) 


r dA "I aQ 

m Jx=xo W 


+ 


ag 

ax 


o 


(3.7) 
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The Eqt3.7) is a first order partial differential equation which is 

non-linear in so far as f~l is a function of discharge 

^dQJx=xo 

Let AB be a river reach for which the stage -discharge 
relationship is assumed to be one-to-one such that Eq(3.7) can be 
applied. 



Fig. 3. 1 A river reach 


Then as discharge Q solely depends on distance x and time t. 


the following relationship holds good: 

dQ = - 3 — dx + ST dt 
da di 


C3.8> 


In order to deternune the constant -discharge propagation lines 
dQ must be zero, hence the following relationship 


df 


i-dQ/d\> 

LdQ/dx) 


dQ 

dA 


(3.9) 


Eq(3.9) is the characteristic equation for the Eq(3.7). 
The discharge (Q) propagation rate is given by 

^ — £13. 

^ dt dA 

If constant celerity c is assumed irrespective of 


Eq(3.7) 


is linear as follows 
1 gQ . ^ ^ 
c gt gx ~ 


0.0 


3 -s- = c = constant 
dA 


(3. 10) 
Q, the 

(3. 11) 

(3. 12) 


or 


gQ gQ 

gt gx 


0.0 


(3. 13) 
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3.2 LEAP-FROG FINITE DIFFERENCE SCHEME: 

The numerical solution of Eq(3.13) can be attempted by a 
finite difference approach assuming a finite difference in the 
CXrt) plane as shown in Fig. 3.2. 



a Leap Frog Scheme 

The partial differentials are approximated by the Leap Frog 
finite difference scheme a^ : 


dQ 

at 

£Q 

dx 


Q"*** - O’?-" 

J 3_ 

2 At 
- Q” 

j+i j-i 

2 Ax 


(3.14) 


If the resulting finite difference equation is 
approximation of differential Eq(3. 13) and if 
con'i/erges with that of Eq(3. 13) when Ax and At tend 


LQ = 


dQ 

at 


+ c 


ao 

ax 


= 0.0 


an acceptable 
its solution 
to zero. 




O’?"" - 
j J_ 

2 At 


+ c 


Q . - Q . 

J+i J-i _ 

2 Ax 


0.0 


(3.15) 


where L, is differential operator based on the Eq(3-14) scheme. 

n 
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On the assumption that these grid function Q(JAx,nAt) can be 
expressed as a Taylor series about point (JAx,nAt) it can be shown 
that the approximation error R is determined by the following 
formula 


R = 1 

R - -g- 


^ 2 2 
1 - r c 




dQ*- 


ax" 


Ax + Higher order terms 

(3.16) 


where r = 


At 

Ax 


The order of approximation by using the scheme of Eg (3. 5) is 2- 

C R = 0(Ax*) 3 

From the Eq(3.15) the following is obtained. 

c At 


— n-i-1 _ n-1 

Q . = Q. 

J J 


Ak 




- Q” 


(3. 17) 


The proposed model , which is based on a Leap-Frog finite 
difference scheme in the MC method, is known as 
Muskingum-Cunge-Leap-Frog (MCLF) model in the rest of the chapter. 


3.3 STABILITY ANALYSIS : 

For the stability analysis in the continuity equation , the 
finite difference form of differential equation can be written of 
the following form. 


_ Q* ^i{/3nAt + C*jAx> 
j 

_n+l i</5<n+l>Al ajAx> 

Q , — u e 

j • 


(3.18a) 


(3. 18b) 
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.= Q* </3<n~i>Ai * o^Ax> 
j 


(3.18c) 


Qr> _ g* ^v</?r»Al 4- 0(Cj4-i>Ax> 
j4-l 


(3.18d) 


qH _ g* ^i</?nAt ■*• OKj-l>Ax> 
j-1 


(3.18e> 


Substituting these set of EqsO.lS) in EqO.lS), and 
re-arranging we have 

gn^gi/JAl _ ^-i^Alj ^ CM gn ^giOlAX_ ^-idAXy ^ 0 g 

J Ax j 

If we define CN = and solving for e'"^^ • 


zvckAx 

e 


ioAx 

e 

"CFT" 


^^i/3At_ g-i/5Alj 


1.0 s 0.0 


The expression for e* ** can be obtained by solving Eq(3.20). 

In the resulting expression substitute a = ia^. Equating the 

real and imaginary terms on both the sides in the expression, it is 

found that for the real term, c* value lies between 0 and 1.0 for 

2 

all values of /3At , which implies that the scheme chosen is stable. 


3.4 Method of Routing: 

The given reach length L, is sub-divided into n number of 
reaches depending on the Ax chosen. In the MCLF method a minimum 
number of three reaches are required. The discharges at all times 
at the next sub reach from the origin are obtained by using the ilC 
method as explained in Section 2.2 in Chapter 2. Then Eq(3.17) is 
used to get the discharges at all reaches with respect to different 


t imes 
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Case 1: 

In this study, a inflow hydrograph t IH) of data set3(Fig. 
3.3) is considered. The rectangular channel has a width of 50 m. 
The given IH is routed over a reach length of 18000 m, on a 

longitudinal bed slope of 0.0005. The lianning^'s roughness 

coefficient was chosen as n = 0.03. A Courant number of 0.8S is 
used for the MCLF and MC methods. A value of Ax = 4500 was used for 
the MCLF and MC methods. The given IH is routed by the MCLF and the 
MC methods. The outflow hydrographs obtained by the MCLF and the MC 
methods in this study are plotted in Fig. 3.3 along with the 

corresponding IH. Fig. 3.4 is an enlarged plot to show the 

attenuation of peak discharge by the MCLF and the MC methods. Some 
of the salient routing data are tabulated in Table 3.1. 

From Figs. 3.3 and 3.4 and table 3.1 it is found that: 

1. The attenuation of peak discharge is very little in MCLF 

I 

method. 

2. The wave translation of the peak discharge is of the 
same order as that of the MC method. 

Case 2: 

An IH of data set 4 CFig. 3.5) with a a steep rise in the 
initial flows was considered next. The physical characteristics of 
the channel are same as that of case 1. A Courant number of O.SS is 
used for the MCLF and the MC methods. A value of Ax = 4500 was used 
for the MCLF and the MC methods. The IH of data set 4 is routed by 
the MCLF and the MC methods. The outflow hydrographs obtained by 

'v. 

the MCLF and the MC methods are shown in Fig. 3.5 along with the 
corresponding IH. Fig. 3.6 shows the attenuation of the peak 
discharge in MCLF and MC methods. Some of the salient routing data 



are presented 

Table 3.1 

in Table 3.2. 

SOME SALIENT 

ROUTING DATA FOR 
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DATA SET 3 

Time 

(hrs) 

ord. of 

Input hydrograph 
m^'/s 

Ord. of output 

MC 

hydrograph (m*/s) 
MCLF 

0.0 

10.0 

10.0 

10.00 

2. 0 

13. 0 



2.28 


10.68 

10.16 

4. 0 

SO. 0 



4.56 


21.77 

19.32 

6. 0 

107. 0 



6.08 


45.77 

40.14 

3. 0 

147. 0 



8.36 


105. 10 

106.90 

9.0 

150.0 


• 

iO. o 

146. 0 



10.64 


144.24 

147.35 

11.4068 


146.259* 

149.9972 

12 . O ■ 

105. 0 



12. 16 


142.93 

147.41 

14. 0 

59.0 



14.44 


102.51 

104.45 

le. O 

33. 0 



16.72 


55.88 

52.57 

13.0 

17. 0 



18.25 


35.91 

33.86 

20. 0 

10. 0 



20.53 


17.21 

16.53 

22.0 

10. 0 



22.05 


11.47 

10.56 

24.0 

10. 0 



24.33 


10.03 

10.00 

25.85 


10.00 

10.00 

. 23.0 

10.0 


. 

28.13 


10.00 

10.00 

* represents peak values. 



From the Figs 3.5 and 3.6 and Table 3.2 


the same 


observations of case 1 are noted. 


Table 3.2 

SOME SALIENT ROUTING DATA 

FOR DATA SET 4 

Time 

Chrs) 

Ord. of 

Input hydrograph 
m^/s 

Ord. of 

MC 

output hydrograph (m^. 

MCLF 

0.0 

10.0 

10.0 

10.00 

2. 0 
2.28 

B2. 0 

11.68 

10.77 

4. 0 
4.56 

141.0 

48.93 

52.34 

6. 0 
6.08 

144. 0 

122.74 

^ 125.52 ^ 

7.48 


142.96 

149.965 

e. 0 

8.36 

114. 0 

139.78 

143.96 

10. 0 
10.64 

86. 0 

121.44 

110.08 

12.0 

12. 16 

61.0 

88.44 

90.51 

14.0 

14.44 

41.0 

62.54 

60.91 

le.o 

lb,72 

25. O 

39.88 

37.31 

IB. O 
18.25 

14. 0 

24.91 

24.06 

20.0 

20.53 

10.0 

14.21 

13.73 

22.0 

22.05 

10.0 

11.47 

10.08 

24. 0 
24.33 

10.0 

10.03 

10.00 

25.85 

4 

10.00 

10.00 

20. 0 

2S. 0 

28. 13 

10.0 

10.0 

10.00 

10.00 


* represents peak values. 



Discharge im /s) 



Fi^. 3>6 Attenuation of peak discharge in 


WC and MCLF methods - Steep rise IH 
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3. 5 DISCUSSION: 

From this study as above it has been observed that the MCLF 

method creates very /little attenuation of peak discharge and the 

translation of peak discharge is of the same order as the MC 

method. In the MC ^ method, the error is introduced into the 
/ 

analytical solution of the continuity equation by the assumed first 
order finite difference scheme. In the actual flood wave there is a 
attenuation in the peak discharge due to the frictional effects. 
The numerical diffusion caused due to the finite difference 
operation fairly matches with that of the actual flood wave. In the 
MCLF method, in which there is a second order accurate finite 
difference scheme, the numerical diffusion is eliminated 
considerably and the solution approaches the analytical solution of 
the kinematic wave equation. Thus the attenuation of peak discharge 
in the MCLF method is minimal. This is in agreement with the 
CungeC6] statement wharvaver <a closar approximation of tho 
di/foi'&ntial equation is obtained, this sliminatss ths wawa 
attsnxuxtion Cor amplifications as ths cass may b&^ . Thus, the 

higher order schemes are not desirable in the MC method and no 
improvement of the MC method can be achieved by this way. 



CHAPTER-4. 


MODIFICATION OF MC METHOD THROUGH SEGREGATION OF INPUT HVDROGRAI 

4.0 INTRODUCTION: 

To improve the Muskingum-Cunge (MOmethod without losing its 
simplicity and to simulate the flows which arc considerably away' 
from the range of peak discharges with more accuracy, a segregation 
technique of the inflow hydrograph for use in the MC method is 
proposed. The number of segregations depends on the range of the 
discharges and shape of the inflow hydrograph. The results from the 
proposed method are compared with those of the Preissman and MC 
methods. 


4. 1 EARLIER STUDIES: 

Cunge 161 obtained 'the following expressions for the routing 
parameters K and 6 as function of the physical character ist ics of 
the channel and inflow hydrograph as 


K = 


Ax 


Q 

and e = 0.5 (1 - a x r- > 

B <c> s Ax 
o 

The outflow hydrograph is obtained 
hydrograph by using these parameters over 


(4.1) 

(4.2) 

by routing the 
a reach length 


inflow 
L. The 


procedure of routing is explained in Section 2.3 


of Chapter 2 
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In certain cases, such as ti) in Fig. 4.1 , in which there is a vast 
difference between the initial flows at which the flood starts to 
build up and the peak discharge it attains and Cii) in Fig. 4.2 , 

when there is a rise after the flood started receding , it is 
felt that the determination of K and © by considering only the peak 
discharge as taken by CungeC63r may not simulate the actual 
conditions in the downstream adequately. To overcome this 

shortcoming, PonceClSl has developed the variable parameter model 
in which the routing parameters K and © are calculated for each 
grid cell. This procedure requires data which outwits the simplicity 
of the MC method. 

Kar megamC 123 , in his stratified approach , has divided the 
inflow hydrograph into several strata as shown in Fig. 4.3. The 
relevant routing parameters K and © for each strata are calculated 
the help of past inflow and outflow data through an 
optimization technique. 



It is felt that this approach of KarmegamC123 is physically 
unconvincing and is purely an empirical approach. Further, by using 
the optimization technique, the simplicity of MC method is 
sacrificed and the routing becomes more complex than what is 
^larran led . 



Time (hr s> 


Fig, 4,1 Inflow hydrograph 




4.2 SEGREGATION MODEL 
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To overcome these shortcomings of PonceClSI and Karmegain[12!I r 
the following Segregation model is proposed. The inflow 
hydrographCIH) is divided into n number of segments depending upon 
the range of values. It is proposed that IH be divided into three 
or five segments by chosing the discharge(s) which is(/are) less 
than peak discharge as a basis for segregation. 

Consider Fig 4.4ar where the inflow hydrograph is divided into 
three segments. The discharge Q which is 0.5 times the peak 

50 

discharge , is taken as the basis for segregation. The resulting 
hydrograph now consists of three segmented hydrographs namely OBE, 
EBPCF and FCD. It is noted that the continuity between the 
hydrographs is not lost. Now it is assumed that c ^ F(t) and is a 
constant irrespective of the discharge. This ensures the continuity 
equation to be linear. It may be noted that this is the basic 
assumption of the MC method. 



With the above assumption, the routing parameter K is the same 
for all the three segregated hydrographs OBE, EBPCF and FCD , i.e. 




(4.3) 
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HowBver, the parameter & for each segreQated unit is proposed to be 
calculated as: 

Q 

For segment OBE 0^= 0.5 Cl - g-- - ) (4.4a) 

o 

Q 

For segment EBPCF 0^= 0.5 Cl - g-^ — ) C4.4b) 

o 

Q 

For segment FCD 0^= 0.5 Cl - g-^ — ) C4.4c) 

o 

Where Q = 0.5 Q . 

50 p 

The IH is routed by MC method Ca procedure in Section 2.3) 
after choosing the appropriate value of © depending upon the class 
to which the discharges belong , after they are segregated. 

Fig 4.4b is a complex hydrograph with two peaks. Even in this 
case, the IH is divided by considering Q as the basis for 

50 

segregation and the values of © are calculated by Eqs. C4.4-a,b,c) . 

The proposed method of segregation and routing with the MC 
method is designated as Segregated -Muskingum-CungeCSMC) model in 
the rest of this chapter. 

4.3 Numerical experimentation: 

The efficacy of the proposed SMC method is compared with the 
implicit Preissmann scheme and the MC method through numerical 
experimentation. Five inflow hydrographs Cset3,set4,set6,set7 and 
setSl were considered. For each IH, ca — is used in the 

Preissmann scheme and a Courant number of 0.85 is used in MC and 


SMC methods. 



Discharge* 
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Casel: 
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First, the IH of data set 3 was considered. This IH was routed 
in a rectan 9 ular channel of 50 m wide for a reach length of 18000 m 
on a bed slope of 0.0005. The Manning'^s roughness coefficient was 
chosen as n = 0.03 . The routing of IH was done by the SMC,MC and 
the Preissmann methods. A value of Ax = 6000 m was used in both MC 
and SMC methods. The IH of data set 3, and the outflow hydrographs 
obtained through these three models in the study are shown in Figs 
4.5 and 4.6. An enlarged plot of values predicted upto Q along 

SO 

with the corresponding IH is shown in Fig. 4.6 . Some of the 

salient routing data are presented in Table 4.1. 

From the Figs. 4.5 and 4.6 and Table 4.1, it is seen that s 

1) The flows predicted by the SMC method are slightly 

better than those of MC method up to Q . Beyond Q , both 

so so 

the SMC and MC methods give essentially identical results. 

21 The attenuation of the peak in the SMC method is of the 
same order as that of the MC method. 

31 The wave translation of the peak discharge for the SMC method 
is of the same order as that of the MC method. 

2 ; 

The second IH used in the study was of data set 4 (Fig. 4.71. 
The physical characteristics of the channel are the same as in 
Case 1. A value of Ax = 6000 m is used for the SMC and MC methods. 
The IH of data set 4 was routed using the SMC, MC and the 
Preissmann methods. The outflow hydrographs obtained through these 
three models are plotted in Figs. 4.7 and 4.S. Fig. 4.8 is an 
enlarged plot of the flows upto Q along with the corresponding 
IH. Some of the salient routing data are tabulated in Table 4.2. 
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TableA. 1 SOME SALIENT ROUTING DATA FOR DATA SET 3 


Ord. of inflow Ord. of outflow hydrograph in m*/s 


Time 

hydrograph 

MC 

SMC 

Preissmann 

(hrs) , 





0.0 

10.00 

10.00 

10.00 

10.00 

2. 0 
2.04 

18. 00 

10.33 

10.09 

9.97 

4. 0 
4.08 

50. 00 

17.58 

16.59 

11.00 

6. 0 

6. 12 

107. 00 

46.78 

45.07 

38.74 

8. 0 
8.16 

147. 00 

99.87 

99.67 

109. Ot 

iO. 0 

146. 00 



144.21 

10.20 


140.60 

140.59 


11.0 

12Q. 00 


A 

144.44 

11.56 

12.0 

105. 00 

146.015* 

146.014 

135. 62 

12.24 


142.24 

142.24 


14.0 

59.00 



99.65 

14.29 

16.0 

33. 00 

106.17 

106. 16 

65. 04 

16.33 

18.0 

17.00 

62.48 

61.42 

41.89 

18.37 

20.0 

10. 00 

34.65 

34.00 

26. 71 

20.41 

22.0 

10. 00 

« 17.91 

17.41 

17. 12 

22.45 


10.86 

10.41 


23.81 

24.0 

10. 00 

10. 11 

10.02 

12. 12 

25.85 


10.00 

10.00 


26.0 

10.00 



10.51 

27,8^ 


10.00 

10.00 


28.0 

28.13 

10.00 

9.96 

9.99 

to. 10 


« reprasants peak value 
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Tablet. 

2 SOME 

SALIENT 

ROUTING DATA 

FOR DATA SET 4 


Ord. of inflow 

Ord. of outflow hydrograph in m^/s 

Time 

hydrograph 

MC 

SMC 

Preissmann 

Chrsi 

Cm®/s) 




0.0 

10.00 

10.00 

10.00 

10.00 

2. O 

52. 00 



9.78 

2.04 


11.62 

10.43 


4. O 

141.00 



38. 14 

4.08 


47.51 

43.92 


6. O 

144. 00 



133. 13 

6. 12 


122. 75 

122.52 


7.0 

130. OO 



140.85* 

7.48 


142.96 

142.61* 


e, o 

' 114.00 



136. 66 

8. 16 


141.30 

141.29 


iO.O 

86.00 



111.09 

10.20 


117.14 

117. 14 


i 2 .o 

61.00 



86. 04 

12.24 


88.46 

88.44 


i4. O 

41.00 



64. iO 

14.28 


62.78 

62.54 


16.32 


41.78 

41.68 


18.0 

14. OO 



32.33 

18.37 


25.72 

25.44 


20.0 

10.00 



22.07 

20.41 


14.79 

14.40 


21.77 


11.20 

10.82 


22.0 

10. 00 



15. 06 

23.81 


10.06 

10.00 


24. 0 

10.00 



11.48 

25.85 


10.02 

10.00 


26. O 

10. OO 



10. 35 

27.89 


10.00 

10.00 


28.0 

10.00 



10.07 


» represents the peak discharges. 


From the Figs. 4.7 and 4.8 and Table 4.2, it is seen that: 

The relative performance of the SliC, MC and the 

Preissmann methods are the same as in Case 1. The 
, three conclusions drawn in case 1 are exactly applicable 




here also 


Discharge (m”/s> Discharge (in*/s) 
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Fig. 4,5 Comparison of SMC model with the MC 
and Preissmann models Normal IH 

Xi 



Preissmann models 


Normal IH 


Discharge (■ /s) Discharge (m®/ 5 > 



Fig. 4.7 Comparison of SMC model with the MC 

and Preissmann models Steep rise IH 



Fig. 4.8 Flows upto Q in the SMC, MC and the 

90 

Preissmann models Steep rise IH 
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Case 3: 

Consider IH of data set 6 as shown in Fig. 4. 9. The salient 
feature of this IH is the presence of a second peak. The magnitude 
of second peak discharge is greater than G but less than the 

50 

fi^i'lisst peak . The rise in IH towards the second peak occurs 
before the commencement of the Q in the depletion curve. This IH 

SO 

was routed by the SMC, MC and the Preissmann methods in a 
rectangular channel of 50 m wide for a reach length of 30000 m. The 
Longitudinal bed slope of channel is 0.0005 and n = 0.03 . A value 
of Ax = 7500 m was used for the SMC and MC methods. The IH of data 
set 6, and outflow hydrographs obtained by these three methods are 
shown in Figs 4.9 and 4.10. Fig. 4.10 is an enlarged plot of the 
initial flows predicted upto G , along with the corresponding IH. 

SO 

Some of the salient routing data are presented in Table 4.3. 

From Figs. 4.9 and 4.10 and Table 4.3, it is seen that: 


U The 

flows predicted 

by 

the 

SMC are 

slightly 

better 

than 

that of MC model 

upto 

Q . 
so 

Beyond Q 

50 

both SMC 

and MC 


methods give identical results. 

2) The attenuation of peak in the SMC method is of the same 
order as that of MC method. 

3) The wave translation of peak discharge for SMC method is ' 
of the same order as that of MC method. 

4) The predicted second peak in the SMC method is slightly', 
highart sss 13 X ) than that by the Preissmann method.' 

I 

However, it is of the same order as that obtained by the MC 


method 
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Table 4.3 SOME SALIENT ROUTING DATA FOR DATA SET 6 


Time 

(hrs) 

Ord. of inflow 

hydrograph 

Cm /s) 

Ord . of 

MC 

Outflow hydrograph in m^/s 
SMC Pre issmann 

0.0 

10.00 

10.00 

10.00 

10.00 

2. 0 

65.28 



10.63 

2. A3 


10.56 

10.03 

• 

4. 0 

154.72 



6. 07 

4.05 


22.25 

17.13 


6. 0 

168.25 



95. 77 

6.49 


98.59 • 

97.01 


8.0 

154. 72 



159.55* 

8.11 


147.99^ 

147.77^ 


9.73 


162.73 

162.71 


10. 0 

130.00 



154. SI 

10.54 


161.15 

161.14 


12. O 

A 40. 60 



136. 14 

12. 17 


148.93 

148.93 


14. 0 

85.28 



128. 38 

14.60 


128.83 

128.83 


16.0 

43. 48 



109. 44 

16.22 


121.16 

121. 12 


17.84 


86.61 

85.45 


18. 0 

17.92 



74.78 

20.0 

10.00 



. 48. 75 

20.28 


38.82 

37.23 

' . 

21.90 


, 20.41 

18.81 


22.0 

10.0 



31.22 

24.0 

10.0 



20.04 

24.33 


10.95 

10.38 


25.96 


10.13 

10.02 


.26.0 

10.0 



14.29 

28.0 

10.0 

' 


11.73 

28.39 


10.04 

10.00 



» represents the peak discharges. 











Case ,4 
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The IH of data set 7 is shown in Fig. 4.11. In this IH , there 


is a second peak of magnitude > Q and less than the first peak 

50 

which occurs after the Q in the depletion curve has been reached. 

50 

This IH was routed in the same rectangular channel as in Case 3 and 
the results are shown in Figs. 4.11 and 4.12. Some of the salient 
routing data are presented in Table 4.4. 

Table 4.4 SOME SALIENT ROUTINS DATA FOR DATA SET 7 

Ord. of inflow Ord. of outflow hydrograph in m^/s 


Time 

thrs) 

hydrograph 

Cm*/s) 

MC 

SMC 

Preissmann 

0.0 

10.00 

10.00 

10.00 

10.00 

2. 0 

65. 28 



10. 63 

2.43 


10.56 

10.03 


4. 0 

154. 72 



6. 07 

4-05 


22.25 

17.13 


6. 0 

168. 25 



95.77 

6.49 


98.59 

97.01 


8. 0 

154.72 



159.54!t.> ' i 

8.11 


147.99^ 

147.77^ 

\ 

9.73 


162.73 

162.71 

154.61UB>ii 

iO.O 

130. 0 



10.54 


161.15 

161. 14 

135. 50 

iZ.O 

98. 36 



12.16 


148.88 

148.88 


14.0 

65. 28 



112. 61 

14.60 


117.08 

117.08 


16.0 

36.47' 



86.34 

16.22 


91.36 

91.32 


17.84 


65.78 

65. 17 


13.0 

70. 90 



62.26 

20.0 

65. 75 



54. 24 

20.28 


53.46 

49.41 


21.90 


71.54 

73.03 

71.57 

22.0 

30.60 



24.0 

10.00 



59.75 

24.33 


56.72 

57.84 


25.96 
26. 0 

10.00 

30.97 ■ 

29.61 

40.05 

28.0 

10.00 



24. 54 

28.39 


12.14 

10.91 





represent peak value 
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From Figs- 4.11 and 4.12 and Table 4.4, it is seen that: 

The first three conclusions arrived in Case 1 ( and also ii 
Case 2 and Case 3) also hold good here. Further, the seconc 
peak is predicted well by all the three methods viz., SMC , MC 
and Preissmann. 


Case 5: 

In IH of data set 8 (Fig. 4.13) , is essentially same as thai 
of the IH of data 7 (Case 4^ , except that the magnitude of th€ 
second peak is less than the Q . The channel data are the same as 

50 

in Case 4. The results of the routing of this IH are shown in Fig. 
4.13, 4.14 and Table 4.5. 


Table 4.5 SOME SALIENT ROUTING DATA FOR DATA SET 8 


Ord. of inflow Ord. of outflow hydrograph in m“/s 


Time 

(hrs) 

hydrograph 

(m®/s) 

MC 

SMC 

Preissmann 

0.0 

10.00 

10.00 

10.00 

10.00 

2. 0 

65 . 28 



10.63 

2.43 


10.56 

10.03 


4 . 0 

154 . 72 

. 


6.07 

4.05 


22.25 

17.13 


‘ 6. 0 

168.25 



95.77 

6.49 


98.59 

97.01 

* 

8.0 

154 . 72 



159.54 

8.11 


147.99^ 

147.77. 


9.73 


162.73 

162.71 


10.0 

130 . 0 



154 . 61 

10.54 


161.15 

161. 14 


12.0 

98 . 36 



135 . '50 

12. 16 


148.88 

143.88 


14.0 

65.28 



112 . 61 

14.60 


117.08 

117.08 


16.0 

36.47 



36.34 

16.22 


91.36 

91.32 


17.84 


65.78 

65. 17 

• 

13.0 

65 . 90 



62.25 










Discharge (.m /s-» 
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20.0 

20.28 

55. 75 

48.64 

44.07 

49. 83 

21.90 

22.0 

20. 60 

62.91 

65.53 

62. 44 

24. O 
24.33 

iO. 00 

42.95 

43.06 

50. 20 

25.96 

26. O 

iO. 00 

21.72 

19.97 

33. 08 

2S. 0 
28.39 

iO.O 

10.86 

10.24 

20. 56 


* represents peak value. 

From Figs. 4.13 and 4.14 and Table 4.4, it is seen that; 


The first three conclusions arrived in Case 1 (and also in 

Case 2 and Case 31 also hold good here. Further, 

1. The predicted magnitude of second peak by the SMC method is 

slightly more in comparison to the Preissmann method and 
the MC methods. 

2. The second peak predicted by the MC method lies between 
those of the SMC and Preissmann methods. However, the 
percentage over prediction by the SMC relative to the 
Preissmann is less than 5 X , and relative to the MC 
method is than 4 X . 

4.4 CONCLUSION: 

A method of Segregating the IH and using MC method, called 
as SMC is proposed. 

1. This method predicts in all type of hydrographs , the 
flows upto Q more accurately than the MC method. 

!50 

2. yhen the IH is a single peaked normal hydrograph , the 
flows beyond G are predicted to the same accuracy as 
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the MC method. 

3. When the IH is a complex hydrograph with double peaks , 
SMC method gives the outflow hyrograph values comparable 
to that one predicted by the MC method with errors less 
than S y. . 

However, when there is a second peak, , both the SMC and the 
MC methods over predict slightly ( S 10 Z ) magnitude of the second 
peak in comparison to the results of the Preissmann method. 



CHAPTER-5 


LAPLACE TRANSFORM METHOD 

5.0 INTRODUCTION; 

An approach for routing the •flood by a lumped model is 
developed by making use of the solution to a step input function of 
Muskingum equation with the boundary condition QiO*) = 0.0. The 

parameters are estimated as in the Musk ing urn -Cunge(MC) method, by 
relating them to the physical and morphological characteristics of 
the channel. The relative accuracies of the proposed model, is 
compared with the MC method and as well as with the implicit 
Preissmann scheme. 




with the initial conditions 
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1(0) = 0(0) = I. 


and obtained the solution as 
Q = I 




(l-d)“^ exp C 


K(i-e) ■'j 


(5.4) 


Gill C103 obtained the solution to the step input by applying 
the Laplace transform to Eq(5.3) to get 

Q + K(l-©)Cp 0 - Q(0'^)3 = f- KeCpf - 1(0*)] (5.5) 

in which Q and I denote the Laplace transform of Q(t) and I(t) 
respect ively. Here, p is a variable used in transforming F(t) in the 
time domain into F(p) in the Laplace transform domain and 0(0*) and 
l(0*)are the initial values. 0 and I approach the origin (t = 0) 

from the positive side of the axis. The specification of the 

« 

initial conditions in this manner is crucially important for 
situations in which a discontinuity exists at t = O ; Eq(5.5) is 

then solved by GillClO] to give 

n _ i - K eCpI - KO*)] + K(l-©) 0(0*) ,, 

0 i-«: kptl-d) 

Inversion of Eq.(5.6) leads to solution of Q(t). 

For a step input at t *= O, mathematically 
I(t) = O , t < 0 

I(t) = I , t > O (5.7) 

o 

There is a discontinuity in the value of I(t) at t =0, and it 

should be noted that 1(0*) = I . The discontinuity at the upstream 

o 

end is progressively diffused and at the downstream end the flow 
can be taken equal to zero at t = 0 so 'hat 0(0*) = O. 
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ftlso, the Leplace transform of unit step of I at t = 0 is given dy 

I ® 

i= -f- 

Substituting I r and QCO*) values in EqC5.6}, we obtain 

I 

- o 

° “ pci+ kp( 1-©)3 ts.s) 


Eq(5.8) is solved directly to give 


Q = I. 


iu^4 


(5.9) 


The Fig. 5.0 shows the definition sketch of a step input and 

the obtained output for the boundary conditions of ICO"^) = I and 

o 

QCO"") = 0. 

PerumalClTl, Dooge et. al C9!l have taken the initial condition 
for Eq. (5.7) as 

Q<0*) = - ^ 15.10) 

and obtained the solution in the same form as Eq.(5.4),viz 


Q = I 




( 1 - 0 ) 


-X 


exp C - 


K(l-0) 


') 


(5.4) 


The solution proposed by Sill ClOl with Q(0'*')=0 is considered 
superior, even though the boundary condition at t = 0 could be 

disputed C9,173. 

The Bill solution for step input as 



exp C - 


K(l-e) 


') 


(5.9) 


is used in this Chapter to' develop a model for routing an inflow 


hydrograph ( IH) 
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5.2 ROUTING OF INFLOW HYDROGRAPH (IH) BY STEP INPUT FUNCTIONS: 

The given IH is approximated by a stair case function as in 
Gille et. al 1113. 



Consider an I'^tt) in Fig. 5.1 an approximation of I(t). 

The IH at any time can be split and written as 

I^(t)= I(0> + CI(l) - IC0)3 u(t-At) + CI(2) - 1(1)3 u(t- 2At) + 

N 

= 1(0) + .t CI(j) - I(j-1)3 u(t- jAt) (5.11) 

where u(t) is a unit step function. 

The transfer function 

H(I ) = I 
o o 

H CM u(t-nAt)3 = M Cl- exp(-^i~j^3 u(t-nAt) 
where M is an appropriate constant, which is the differences 
between the two adjacent step input functions. 

Since the Muskingum method assumes linear relationship between 
the storage , input and output hydrographs, the principle of 
superposition theorem can be used to get: 
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= ICO) + CI(1)-IC0)I Cl - exp C- uCl-At) 

+ CIC2)-IC1)3 Cl - exp C uCt- 2At) 

+ CIC3)-IC2)1 Cl-expC- uCt- 3At) + 

The output hydrograph at any time t is calculated as 

Q(t) = ICO) +.|^CICj) - ICj-l)l Cl - expc- 3 uCt-jAt). 

C5.13) 

Eq. 5-13 is the output for the given input stair case function 
I‘'Ct). Knowing the MC parameters K and 0, for any given IH the 
output QCt) can be explicitly determined by using Eq.CS. 13). 

The proposed model, which is routed with series of step inputs 
on the basis of Muskingum equation is designated as 
Musk ingum-Lap lace -Transform CMLT) model in rest of the Chapter. 

5.3 ESTIMATION OF PARAMETERS K AND 0s 

5.3. 1. When the input hydrograph and physical characteristics of 
channel are known:* 

The rectangular channel has a breadth B 'units and a 

longitudinal bed slope of s . The Manning's roughness coefficient n 

o 

. N 

is known. The IH has a peak discharge Q . The IH was routed by MLT 

p 

method through a reach length of L units. In the MLT method, to 
ensure the continuity equation to be linear, it is assumed that 
c vt FCt) and is a constant irrespective of the discharge. To route 
an IH by MLT method the MC coefficients K and 0 are proposed to be 


used. 



Thus by MC method, (Eqs 2.30 and 2.31) 
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~ ^c> ' where K has units of time. 

G 

and 0 = 0.5 Cl - g_g-E_^) and 0 < © < 0.5. 

o 

The value of wave celerity c, could be further assumed by 
either of the following two methods: 

1. Assume a celerity of c value and then iterating as per 
Raudikivi's algorithm given in Section 2.3.2(Chapter2, pagel3) . 

2. The wave celerity c, can be calculated by the 

Kleitz -Seddon's lawCEg. C2.9) 3 as 


c 


C 


dQ 

dA 


) 


X 



C-^) 

dy^ 


X 


dQ 

where = slope of the stage-discharge curve at any 

section x. 


In the absence of the stage -discharge relationship, Manning''s 

dQ 

formula may be used to estimate (.-r—) and hence c. 

dy x 


5.3.2 When the Input and. Output hydrographs are available: 

Uhen one or more sets of inflow and outflow hydrographs are 
available for the catchment, by the method of moments the different 
sets of parameters of K and 0 are obtained for different inflow and 
outflow hygrograph sets Cas explained in Section 2.2.2). 


Thus 


K = M'^CQ) 

X 


M°CI) 

1 


0 = 0.5 


{ M CQ) - M (I) 1 


where, M^CQ) and M°CI) are the first moments of the outflow and 

inflow hydrographs about the origin respectively. 

M CO) and N Cl) are the second momer’.:' 

2 2 


of outflow and inflow 
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hydrographs about their center of gravity respectively. 

From these sets of K and 0 , the appropriate K and 0 are 
chosen for use in the MLT method depending on the range of IH 
values. 

5.4 PROPOSED AL60RITHM OF MUSK INGUM -LAP LACE -TRANSFORM (MLT) : 

The following procedure is adopted for routing a 
flood by MLT method. 

Given ; Rectangular channel of breadth = B , Slope = s 

o 

Manning-'s roughness coefficient = n , Reach length = L 
IH with ordinates at time interval DTG. 


N 

1. Choose a value of DT at which the IH of period DTG, can be 
represented by a staircase approximation function. 

2. Obtain the values of the IH at an interval of chosen DT with a 
a second order interpolation procedure. 

3. From the interpolated IH values , obtain the peak discharge Q^. 

4. Calculate the normal depth y for the peak discharge Q . 

n P 

5. Calculate the wave celerity c for the rectangular channel by 
Kleitz-Seddon's law as: 


where ft is the exponent value for the discharge area relationship 
based on the resistance formula- For Manning^'s formula ft = 1.66. 

6. Calculate the routing parameters K and 0 as : 


K = 


0 = 0.5 (1- 


B s <c>L 
o 


7. The approximated is now routed by the using routing paramo 'ers K 
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and ©, to get the outflow hydrograph ordinate at any time t is 
obtained as : 

QttJ=ICO) + ^r^ C ^l-exp(- j u<t-jAt> 

where uCt) is a unit function. 

5.5 NUMERICAL EXPERIMENTATION: 

The efficacy of the proposed MLT method is compared with the 
implicit Preissmann scheme and the MC method through numerical 
experimentation, using four sets of data Cset2-51. For all sets of 
data, a value of a = 0.66 is used in the Preissmann scheme and a 
Courant number of 0.85 is used in the MC method. The MLT model was 
approximated into a staircase function at a, time step DT = 100 

seconds <DTG/DT = 36). Each IH was routed by the MLT, MC and the 
Preissmann methods. 

Case 1 : 

First, the normal IH of data set 3(Fig5.2) for a rectangular 
channel was considered. The channel has a width of 50 m and a 
longitudinal bed slope' of 0.0005. The Manning's roughness 
coefficient was chosen as n = 0.03. The IH was routed over a reach 
length of 18000 m. A value of Ax = 4500 m was used in the MC 

method. The outflow hydrographs obtained by the MLT, MC and the 
Preissmann methods in this study are plotted in Fig. 5.2 alongwith 
the corresponding IH. Fig. 5.3 is enlarged plot to show the 
attenuation of the peak discharge in these three methods. The 
numerical error is calculated for the peak discharge and its time 
of occurrence as: 

E - E 

y. Error = — g — - X 100 
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where is the Preissmann peak discharge or its time of 
occurrence. For all the cases (case 1-4), these error percentages 
arc tabulated in Table 5.5. Some of the salient routing data for 
the MLT, MC and the Preissmann scheme are presented in Table 
5.1Ca). The peak characteristics obtained by the MLT, MC and the 
Preissmann methods arc tabulated in Table 5.1(b). 

From the Figs. 5.2 and 5.3 and Table 5.1(b) it is seen that: 

1. The attenuation of the peak discharge in the MLT method is 
slightly more (1.23)i) in comparison to the Preissmann scheme. 

2- The attenuation of the peak discharge in the MC method is 

slightly less (1.26X) in comparison to the Preissmann scheme. 

3. The translation of the peak discharge in the MLT method 
is slightly less C6.62X) in comparison to the Preissmann 
scheme. 

4. The translation of the peak discharge in the MC method is 
slightly more(3.69X) in comparison to the Preissmann scheme. 

5. The absolute error in the prediction of the peak 


discharge by the MLT metTiod is of the same order as that 
of the MC method. 


6. The absolute 

error 

in 

the prediction 

of 

t ime 

of 

occurrence for 

the 

peak 

discharge is 

more 

in 

the 


than the MC mett^od. 



Table 5.1 (a) 


SOME SALIENT ROUTING DATA FOR DATA SET 3 


Time 

Chrs) 

Ord.of 

Ordinates of 

outflow 

a 

hydrograph in m /s 

1 n T 1 ow 
(m^/s) 

Pi^eissmann 

MLT 

MC 

0.0 

10.00 

10.00 

10.00 

10.00 

2.0 

12.00 

9.97 

12.82 

10. 60 ^ £. 28 hr 

4.0 

50.00 

11.00 

29.41 

21.77 ® 4. 56 hr 

£>. 0 

107.00 

38.74 

71.79 

45. 77 & 6.08 hr 

8.0 

9.0 

147.00^ 

150.00 

109.01 

120.98 

105. 10 & 8. 36 hr 

10.0 

146.00 

144.21 

142.09 
142.66 % 

10.27hr* 

11.0 

iSQ. 00 

144.44* 


146.26 e 11.406 hr* 

12.0 

105.00 

135.62 

127. 44 

142.93 & 12.17 hr 

14.0 

59.00 

99.66 

86.78 

102.51 & 14.45 hr 

16.0 

33.00 

65.04 

52.72 

68.95 & 15.97 hr 

18.0 

17.00 

41.89 

29.59 

35.92 @ 18.25 hr 

20.0 

10.00 

26.71 

15.96 

22. 14 & 19.77 hr 

22.0 

10.00 

17.12 

11.37 

11.47 ® 22.05 hr 

24.0 

10.00 

12.12 

10.32 

10. 03 & 24. 33 hr 

26.0 

10.00 

10.51 

10.07 

1 

i 

10. 00 & 25. 85 ?%r i 

28.0 

10.00 

10.10 

10.01 

! 

9. 99 & 28. 13 hr 


» represents peak discharge values. 


Table S.lCb) PEAK CHARACTERISTICS FOR DATA SET 3 


Preissmann MLT MC 


Peak discharge Attenuation . 5.56 7.34 3.74 

tin m^/s) 

Time to peak Cin hrs) 11.0 10.27 11-406 


X)ischdr9C^Rk Discharge (.» / 


iS0 



Fig. 5.2 CompariBon of MLT model with the MC 

and the Preissmann models - Normal IH 



Fig- 5.3 Attenuation of peak discharge in MLT, 

MC and the Preissmann models - Normal IH 
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Case 2: 

A steep rise IHCFig. 5.4) of data set 4 for a rectangular 
channel was considered next. The physical character ist ics of the 

channel are same as in case 1. a value of Ax = 6000 m is used for 

the MC method. The IH was routed by the MLT, MC and the Preissmann 
methods. The outflow hydrographs obtained by these methods are 

platted in Fig. 5.4. Fig. 5.5 is an enlarged plot to show the 

attenuation of peak in the MLT, MC and the Preissmann methods. Some 
of the salient routing data arc tabulated in Table 5.2(a). The peak 
characteristics are tabulated in Table 5.2(b). 

From the Figs. 5.4 and 5.6 and Table 5.2(b) it is seen that; 

1. The attenuation of the peak discharge in the MLT method 
is slightly more(3.33>C) in comparison to the Preissmann 
method. 

2. The attenuation of the peak discharge in the MC method is 
slightly less(1.277.) in comparison to the Preissmann method. 

3. The translation of the peak discharge in the MLT method 

is slightly less(5.42X) in comparison to the Preissiriann method 

4. The translation of the peak discharge in the MC method 'is 
slightly more(6.85X) in comparison to the Preissmann method. 

5. The absolute error in the prediction of the peak 

by the MLT method is more than that of the MC method 
The absolute error in the prediction of peak discharge by 
the MLT method is more than that of the MC method. 


6 . 



Table 5.2(a) 


SOME SALIENT ROUTING DATA FOR DATA SET 4 
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Time 

(hrs) 

Ord. of 

Ordinates 

of outflow hydrograph in m^/s 

I n f 1 ow 
(m^/s) 

Preissmann MLT 

HC 

0.0 

10.00 

10.00 

10.00 

10.00 

2.0 

52.00 

9.78 

24.58 

il.62 €> 2.04 hr 

4.0 

5.0 

141.0 

150.0 

38. 14 

88.94 

47. 5i 4.08 hr 

6.0 

7.0 

144.0 

133. 13 

140.85* 

134.04 

136.16 © 6.62hr* 

122. 75 ® 6.12 hr 

142.64*© 7.48 hr 

8.0 

114.00 

135.66 

128.00 

141.31 @ 8.16 hr 

10.0 

86.00 

111.09 

104.00 

117. 14 @ 10. 26 hr 

12.0 

61.00 

86.04 

78. 17 

88.46 @ 12.25 hr 

14.0 

41.00 

64.11 

55. 12 

62. 79 @ 14.29 hr 

16.0 

25.00 

46.23 

36.32 

41.38 @ 16.32 hr 

18.0 

14.00 

32.33 

22.22 

25. 72 18.37 hr 

20.0 

10.00 

22.07 

13.57 

17.80 @ 19.73 hr 

22.0 

10.00 

15.06 

10.82 

10. 48 <® 22. 45 hr 

24.0 

10.00 

11.48 

10.19 

10. 06 @ 23. 81 hr 

26.0 

10.00 

10. 35 

10.04 

10.02 e 25.85 hr 

28.0 

10.00 

10.07 

10.00 

10.00 @ 27. 89 hr 

* represents peak 

discharge 

values. 


Table 5. 

2(b) 

PEAK CHARACTERS I ST ICS FOR DATA SET 4 




Preissmann MLT 

MC 


Peak discharge Attenuation 9.15 13.84 7.36 

(in iri^/s) 


Time to peak (in hrs) 


7.0 


6.62 


7.48 



Dxschan 








Case 3: 
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A triangular IHCFig. 5.6) of data set 5 for a rectangular ’ 
channel was considered. This IH of Fig. 5.6 has a peak discharge of 
2000 m /s occurring at 10.0 hrs. The channel has a width of 60 m 
and a longitudinal slope of 0.0015. The Manning'^s roughness 
coefficient n = 0.02 was chosen. The IH was routed by the MLT , MC 
and the Prcissroann methods through a reach length of 30000 m. A 
value of Ax = 7500 m was used for the MC method. The outflow 
hydrographs obtained in the study by the MLT, MC and the Preissmann 
methods are plotted in Fig. 5.6 along with the corresponding IH. 

Fig. 5.7 is an enlarged portion of the plot to show the peak 

attenuation in the three methods. Some of the salient routing data 

/ 

are presented in Table 5.3(a). iThe peak characteristics are shown 
in Table 5.3(b). 

From the Figs 5.6 and 5.7 and Table 5.3(b) it is seen that: 

1. The attenuation of the peak discharge in the MLT method 
is slightly more (4.77X) in comparison to the Preissmann 
method. 

2. The attenuation of the peak discharge in the MC method is 
slightly less (3.12%) in comparison to the Preissmann 
method. 

3. The translation of the peak discharge in the MLT method 

is slightly less(1.69%) in comparison to the Preissmann 
method. f - 

4. The translation of the peak discharge in the MC method is 
slightly more(il.77%) in comparison to the Preissmann method# 

5. The absolute error in the prediction of the peak 

discharge by the MLT method is more than the MC method. f 
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absolute error in the prediction of the peak discharge 
by the MLT method is less than the MC method. 


Table 5. 

3(a) 

SOME SALIENT 

ROUTING DATA FOR 

DATA SET 5 

Time 

(hrs) 

Ord. of 

Ordinates 

of outflow hydrograph in /s 

Inflow 

(m^/s) 

Preissmann MLT 

MC 

0.0 

50.00 

50.00 

50.00 

50.00 

2.0 

80.00 

49.89 

60.61 

50. 00 @ 2. i, hi- 

4.0 

110.00 

51.42 

82.99 

st, ii ^ 4.03 hi 

6.0 

140.00 

78.48 

109.98 

75. 88 @ 5. Q7 

8.0 

170.00 

115.52 

138.80 

107. 37 @ 8.07 hi 

10.0 

200.00* 

151.26 

168.34 

136. 44 @ 10. 01 hi 

12.0 

13.0 

185.00 

151.25 

191.97* 

182.66 

182.81 e 12.78 hr 

165.51 @ 11.95 hi 
197.94 @ 14.53 hr* 

14.0 

170.00 

190.41 

178.77 

f95. 75 @ 14.04 hi 

16.0 

155.00 

177. 47 

168.15 

187. 02 & 15.98 hi 

18.0 

140.00 

163.29 

154.67 

171.22 <§> 18. 08 hi 

20.0 

125.00 

149.15 

140.54 

156. 69 @ 20. 02 hi 

22.0 

110.00 

135.12 

125.81 

142. 1'5 @ 21 .96 hi 

24.0 

95.00 

121.23 

110.91 

126.41 @ 24.06 hi 

26.0 

80.00 

■107.52 

95.95 

111.87 

28.0 

65^00 

94.03 

80.97 

96. 13 @ 28. 09 hi 

30.0 

50.00 

80.85 

65.97 

81.59 <® 30. 03 hi 

32.0 

50.00 

68.03 

56. 14 

67.06 @ 31.97 hi 

34.0 

50.00 

55.94 

52.42 

51.79 €> 34.07 hi 

36.0 

50.00 

51.45 

50.94 

49.98 ® 36.01 hi 


» repressnis peak value. 
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Table 5.3(b) PEAK CHARACTERS I ST ICS FOR DATA SET 5 




Preissmann 

MLT 

MC 

Peak discharge 

Attenuation 

8.03 

17. 19 

2.06 

(in m^/s) 





Time to peak (in 

hrs) 

13.0 

12.78 

14.53 


Case 4: 

A Pearson type III hydrograph (data set 2) which has a stage 
hydrograph as shown in Fig- 5.8 for a rectangular channel was 
considered. The channel has a width of 5000 ft (1500 m) and a 
longitudinal bed slope of 0.000189. The Manning^'s roughness 
coefficient n = 0.03. This hydrograph was routed through reach 
length of 100 miles (160 Km) by MLT, MC and Preissmann methods. The 
outflow stage hydrograph obtained in this study are plotted in 
Fig. 5.8. Fig. 5.9 is an enlarged plot of Fig. 5.8 to show the 
attenuation of depth by the MLT, MC and the Preissmann methods. 
Some of the salient routing data are tabulated in Table 5.4(a)- 
Peak characteristics are tabulated in Table 5.4(b). 
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Table 5.4(a) SOME SALIENT ROUTING DATA FDR SET 2 


Time 

(hrs) 

Initial 

(ft) 

Final stage 

hydrograph 

in ft 

Preissmann 

MLT 

MC 

0.0 

5.08 

5.08 

5.08 

5.03 

6.0 

5.22 

5.00 

5. 10 

5. 08 & 5.07 hr 

12.0 

7.21 

5.00 

5-66 

5. i5 <® i2. 69 hr 

18.0 

12.04 

4.98 

8.09 

24.0 

18.02 

5.00 

12.42 

5. 64 @ il. 77 hr 

30.0 

23.49 

5.51 

17.51 

8. 72 & 25. 39 hr 

36.0 

27.57 

2.86 

22.29 

12. 46 @ 30. 47 hr 

42.0 

29.95 

10.20 

26.04 

16.85 ® 35.55 hr 

48.0 

30.68 

21.62 

28.47 

£3. <S 43.17 

54.0 

30.06 

26.85 

29.51 

26. 42 <§> 48. 25 hr 

60.0 

28.43 

29. 

29.01 

585 e 56 hr 
29.34 

* 28. 67 ® 53. 33 hr- 

30.04 @ 60.94 hr* 

64.5 

66.0 

26.15 

29.4 

29.36 

28.20 

29.78 @ 66.02 hr 

28.04 ® 73.64 hr 

72.0 

23.52 

28.57 

26.36 

78.0 

20.80 

27.06 

24.09 

26.26 ® 78.72 hr 

24.20 ® 83.80 hr 

84.0 

IS. 15 

25. 17 

21.62 

90.0 

15.70 

23.09 

19.12 

20.88 & 91.42 hr 

96.0 

13.51 

20.99 

16.73 

18.69 ® 96.5 hr 

102.0 

11.62 

18.97 

14.54 

15. 62 @ 104. i hr 

108. 0 

10.05 

17.07 

12.60 

13. 80 @ 109. 2 hr 

114.0 

8.77 

15.34 

10.93 

12. 17 @ 114. 3 hr 

120.0 

7.76 

13.79 

9.53 

10.76 @ 119. 3 hr 

126.0 

6-99 

12.41 

8.40 

9. 03 ® 126. 9 hr 
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132.0 

6.41 

11.20 

7.50 

8. 12 


132.0 hr 

138.0 

5.99 

10.14 

6.85 

7. 38 


137. 1 hr 

144.0 

5.70 

9.23 

6.30 

6. 55 


144.75 hr 


* represents peak depth value. 


Table 5.4(b) PEAK CHARACTERS I ST ICS FOR DATA SET 2 

Preissmann MLT MC 

Peak depth Attenuation 1.28 1.095 0.64 

(in ft) 

Time to peak (in hrs) 64.5 56.0 60.94 


From the Figs 5.8 and 5.9 and Table 5.4(b) it is seen that: 

1. The attenuation of the peak depth in the MLT method is 
of the same order as the Preissmann method. 

2. The attenuation of the peak depth in the MC method is 

slightly less (2.187.) in comparison to the Preissmann 
method. 

3. The translation of the peak depth in the MLT method is 
slightly less(13. 177.) in comparison to the Preissmann method. 

4. The translation of the peak depth in the MC method is 

slightly less(5.52%) in comparison to the Preissmann method. 

5. The absolute error in the prediction of the peak 

discharge by the MLT method is less than the MC method. 

6. The absolute error in the prediction of the peak discharge 

by the MLT method is more than the MC method. 

0 
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Fig. 5.8 Comparison of MLT model with the liC 

and the Preissmann models - Log Pearson Type III IH 



Fig. 5.9 Attenuation of peak depth in MLT, 

hC and the Preissmann models - Log Pearson type III IH 
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Out of Curiosity, the solution of Dooge.et.alC93 i.e. Eq.(5.4) 
was used instead of Eq.(5.93 in the MLT method. This method is 
referred as DMLT method. The results of "case 1 are compared by 
routing the IH of case 1 by DMLT [Figs. 5.10, 5.111. 

From the Figs. 5.10 and 5.11 it is seen that: 

1. The output hydrograph values obtained in the DMLT 
method CDooge.et.al'^s solution(Eq.5.4) 3 has negative 
flows in the initial region. 

2. The outflow hydrographs values obtained by DMLT method 

Dooge. et . al "'s solutionCEq.S. 43 do not match at all with 
the outflow values of either the Preissmann method or 
the MLT method. 

As such, it is concluded that Gill''s solutionCEq.5.93 is 
superior and relevant than the Dooge. et . al ^s solutionCEq. 5.43. 

) 

5.6 CONCLUSION; 

Using Gill's solution CEq.5.93 a method of routing the IH is 
proposed, called as MLT method. In the MLT method a stair case 
function is used to represent the given IH. An algorithm for the 
MLT method is given. The efficacy of the MLT method is compared 
with the implicit Preissmann and the MC method through numerical 
cxper imentat ion using four sets of data. 

Table 5.5 gives an abstract of the relative performance of 
the MLT, MC and the Preissmann methods of routing for the four 
inflow hydrographstdata set 2- 4). 



Discharge Cm /s> 
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Table 5.5 








Error 7. 

Error 7. 

Data set# 


Preissmann MLT 

MC 

Far MLT 

For MC 




For peak 

dischargc(m^ 

/s)/depth(ft) 


2 

(Normal ) 

Q 

144.44 

142.66 

146.26 

-1.23 

1.26 

3 

(Steep) 

Q 

140.85 

136. 16 

142.64 

-3.33 

1.27 

5 

(Tr iang. ) 

Q 

191.97 

182.81 

197.94 

-4.77 

3. 12 

4. 

(Pearson) 

H 

29.401 

29.585 

30.042 

0.63 

2. 18 




Time of 

occurence (in hrs) 



2 

(Normal ) 


11.00 

10.27 

11.41 

-6.63 

3.69 

3 

(Steep) 


7.00 

6.62 

7.48 

-5.42 

6.85 

4 

(Tr iang . ) 


13.00 

12.78 

14.53 

-1.69 

11.77 

5 

(Pearson) 


64.50 

56.00 

60.94 

-13.17 

- 5.52 


It is found that 

1. The MLT method as proposed is a viable alternative to 
the MC method. 

2. In general, the MLT method predicts peakflow to the same 
degree of accuracy as the MC oiethod. Howevjzr, while the 
MC method was' found to over predict the peak (with errors 
of order of 3'/) relative to the Preissmanh scheme, the 
MLT method under predicts the peaktwith errors of order 
of 5X) relative to the Preissmann scheme. 

3. The time of occurrence in the peak is also predicted 
substantially to the same degree of accuracy as the MC 
method. 

These results are for a wide range of IH^'s , having normal, 
Etesp, triangular and log psarson type III patterns. The MLT method 
gives coroparable accuraci^^s in a majority of cases to those of the 

jiC 
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Special Features Of The MLT Methods 

In addition to the ^ above performance characteristics the 
following general advantages of the MLT method are worth noting: 

1. In the use of the MC methodr the IH has to be routed in 
through various sub reaches in between whereas in the MLT 
method, the IH is routed in a lumped manner over the full 
reach. 

2. TheMC method has to satisfy Courant stability condition 
whereas, the MLT method does not have any stability 
conditions to be taken care off. 

3. The Programming and Computation efforts are considerably less 
in the MLT method when compared to the Preissmann and the 


MC methods 



CHAPTER-6 


CONCLUSIONS AND RECOMMENDATIONS 

6.1 The basic Muskingum method, the Muskingum-Cunge (liC) method and 
the various improvements of MC method are discussed in this study. 
The forward finite difference scheme is changed into a Leap frog 
finite difference scheme, called as the MCLF method. By segregating 
the IH and using the IH and fusing MC method a new model called SMC 
is proposed. An approach for routing the flood in a lumped manner, 
called as MLT method is proposed. 

6.2 The following conclusions are drawn after comparing the MC 
method with the implicit Preissmann scheme; 

1. The MC method over predicts the peak relative to the 
implicit Preissmann scheme. 

2. The time of occurrence of the peak is also over predicted 
in the MC method relative to the implitit Preissmann 
scheme. 

The MC method, predicts the initial flows more than the 
implicit Preissmann scheme. 


3 . 



6.3 The following conclusion 
experimentation on the MCLF model: 


drawn after 


84 


is 


numerical 


1. The higher order schemes are not desirable in the liC 
method and no improvement can be achieved by this way. 

6.4 On comparing the SMC model with the MC and the Preissmann 
methods, by numerical experimentation on five sets of data the 
following conclusions are drawn: 

1. The BMC model predicts the initial flows better than the 
MC method upto the chosen segregating discharge CQ ) for 

- 50 

all inflow hydrographs. Beyond Q , the SMC predicts the 

50 

outflows to the same accuracy as the MC method. 

2. Thus, with the SMC method, the efficacy of the MC method 
is generally improved. 

6.5 The following significant conclusions are drawn after the 
numerical experimentation and comparison of the MLT model with the 
Preissmann methods: 

1. The MLT model predicts outflow peak values and its time of 
occurrence to the same degree of accuracy as that of the 
MC method. 

2. The MLT model has a distinction of being a lumped model 
and it can be programmed easily- There is no restriction 
of Courant conditions on the time step size as in MC 
method. 

The MLT model is an alternate, viable to the MC method. 


3. 
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RECOMMENDATIONS: 

1. Further studies are needed in the MC method to properly 
estimate the routing parameters K and 9 . 

2. To know the behavior of the SMC method for comp lex 
hydrographsr more experimentation is desirable. 

3. As the MLT method is a viable alternative to the MC method, 
it is suggested that further work be carried out in 
the estimation of routing parameters K and 9 applicable to 


the MLT method. 
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APPENDIX-1 


PRIESSMANN SCHEME 
FOR 

NUMERICAL SOLUTION OF St.VENANT EQUATIONS 

"For the hydrodynamic flows the governing Saint-Venant 
equations constitute a set of non-linear, first-order, hyperbolic 
partial differential equations. For the Saint-Venant equations, 
general analytical solutions ar^ not available except for only a 
few idealised cases. Therefore, they are generally solved by using 
numerical techniques. In this section, one such scheme, originally 
developed by Preissmann(1961) is presented. 


~T 

h 


j 

V 

/ 

f 

f i 

I 



^ . 


/ / /"/■’/v"/ / / ;/ 7 y 

h* B ^ 

1 

(a) 

Elevation 


(b) Cross-Section 



Fig. A.l 

Definition Sketch 


In a 

wide rectangular 

i 

prismatic channel as 

shown in 

F 1 g w ^ . 1 , 

the Saint-Venant 

equations are derived 

using the 


following assumptions. 

1. The flow is one -dimensional in the direction of channel. 

2. The channel is straight and prismatic. 
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3. The pressure distribution is hydrostatic. 

4. The velocity distribution at a channel section is uniform. 
S' The transient -state friction losses may be computed using 

formulas for the steady state conditions. 

6. The bed slope of channel bottom is small such that the 
cosine of the angle made by the averaged bed profile with 
the horizontal can be regarded as unity. 

The Saint -Venant equations based on principles of mass and 
momentum to a differential control volume C2!]. The conservative 
form of these equations are : 


m 

dt d'A 


= 0 


m d , Q 
at ax ^ A 


where 


+ g I) = gA (s -s,) 

o 1 

y<x> 

I = / (y - C> <yCC) dC 

o 


(Al) 


(A2) 


in which 


- width of cross-section 
C = depth integration variable 


+ 

ay 


d\ 

at ^ 

dx 

av ^ 

av . 

dy 

at 


dx 


The non -conservative form of Saint-Venant equations are: 

= 0 ’ CA3) 

and ^+v^ + g^ = 9^%“ ) (A4) 

Where y is the depth of flow; v is the flow velocity, 9 is the 
acceleration due to gravity; s^ is the slope of channel bottom; s^ 
is the slope of energy gradient line; x is the distance measured 
downstream direction, and t is the time of computation. s^ is 

evaluated using Manning's formula for a wide rectangular channel. 

2 2 
n V 

y 

where n is the Manning's roughness coefficient. 



PREISSMANN SCHEME: 


92 


Preissm^nn scheme is basically a four -point weighted implicit 
finite difference scheme, used to transform the non-linear partial 
differential Saint— Venant equations into non-linear algebraic 
equations. Referring to Fig. A. 2 a function f(x,t) is represented 

by its values on a grid with f" denoting the value at the grid 

ih ( J ^ 

and n time interval. Ax^. is the non-dimensional grid spacing 

between the + and j^^point. Non-dimensional isat ion is 

explained in the next section. Then for a non-dimensional time step 

At, the four-point weighted difference approximations for the time 

and space derivatives f, and f for a cell defined by the four 

points (n, j) , Cn+1, j) , tn, j+1) , (n+1, j+1) are: 



and 


f (x,t) 


j+-i j 
2 

. 


(1-a) 



(A71 


where a is a weighing factor varying from 0.5 to 1.0. 

Preissmann scheme is first order accurate in space and 
second -order accurate in time- 


BDUNDARY CONDITIONS: 

In the Preissmann scheme, from N nodes (2N-2) algebraic 
equations are obtained from the Saint -Venant equations. The 
remaining two equations are obtained from boundary conditions. In 
the subcritical flows, only one boundary condition is’ sufficient- 
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The downstream can be evaluated by second order extrapolation. The 
upstream boundary condition is specified by a hydrograph. At the 
downstream end it is extrapolated by a second-order extrapolation 
^= 0.0 

These non-linear equations are solved using 
Newton -Raphson methodC23. 


NON -D I MENS I ON AL I SAT I ON ; 

The Saint -Venant equations are non-d imensional ised using 

channel length L as distance parameter and initial steady uniform 

o 

flow depth Y as the depth scale parameter. 

O 

The initial steady state uniform flow velocity, v is 

o 

taken as the velocity scale parameter and is obtained using 



defined as 


X 



o 


<A10) 

(All) 

(A12) 

(A13) 


Substitution of all these Eqs(A8-A13) in Eqs A3 and A4, it results 



in 
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and 


dy 

+ h G 


di 

° dx 

-1- y — 


= 0.0 


(Ai4> 


_i “ 5v « 3y — 

d — -I- e V — + f ^ - p s 3 (A15> 

° ° dx ° dx ° ° ° ^ 

where a^, b^, c^, d^, e^ , f ^ and are dimensionless coefficients. 

is dimensionless constant and it becomes dimensionless when 


combined with 


= 


2-2 
n V 

- 4/3 


tA16) 


The best set of coefficients after non-dimensional isat ion are: 

(A173 


a =1.0 
o 


b = 
o 


V T 
o o 


c = b 
o o 


d = -=r 

o gT 


e = — 
o gL 


(A18) 

tA193 

(A203 

(A21) 


"o = — 
o 

g = 1.0 


(A223 

(A23) 


and 


Po = 


> 4/3 


(A243 


o \ 

FFiOGnAM: A value of a = 0.6& is used in the Preissmann schema. The 
Preissmann schame is programmed and tested against Moin and cooley 
dataCSl. The error between the Moines data and the outflow values 
obtained by the Preissmann schema are in permissible limits. This 
implicit Preissmann scheme is used as a benchmark • to compare the 
proposed models in this study- 



APPENDIX-^ 


INPUT DATA 

The following data sets were used in this study for the 
numerical experimentation. / 

Set # 1. (Cooley and Moin,1976). 

A rectangular channel 20ft (6.1m) wide and 2 mi les(3. 2kni) long 

carrying a steady uniform discharge of 833.49 cfs(23.34 mVs) at 

6ft (1.83m) depth is subjected to an upstream flood - wave with a 

peak of 2000 cfs (56m^/5) increasing linearly in a period of 20 

min. This upstream flow decreases linearly from its peak of 200Ocfs 

to 833. 49cfs(57m^/s to 23.34 m^/s) in 40 min. Additional properties 

given by Veismann .et.al arc s = 0.0015 ft/ft and n=0.02. The 

o 

problem is to follow the translation of the flood wave down the 
stream. 

Set # 2. (Cooley and Moin, 1976) 

A hydrograph described by the Pearson Type III distribution is 
routed down a wide C5000 ft (1500 m) was used herein!- Open channel 
having length of 100 miles(160 km), s^ = 1/5280 ft/ft, and n = 0.03 

The initial depth and discharge are 5 ft (1.5m) and 
49766 cfs (1390m^/s) respectively, and all subsequent discharges 
are given by 




r 

t 

Cl/t^' i)3 Ci-a/T)! 

Q(t)=Q 

o 

l+(p-l) 

T 

e 
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in which p- - 20: t = the time at which the maximum 

occurs; and the initial upstream discharge. 

Set # 3 (Raudikivi, 1979) 

The hydrograph at the upstream end of a river reach(x=0) 
is as follows: 


t hours; 

0 

1 

2 

3 4 

5 

6 

7 8 

9 

Q m^/s : 

10 

12 

18 

28.5 50 

78 

107 

134.5 147 

150 

t hours; 

10 

11 

12 

13 

14 

15 

/ 

16 

17 

Q m^/s 

14£> 

129 

105 

78 

59 

45 

33 

24 

t hours: 

18 

19 

20 

21 

22 

23 

24 

25 

Q m^/s : 

17 

12 

10 

10 

10 

10 

10 

10 

t hours: 

26 

27 

28 






Q m^/s : 

10 

10 

10. 






Length of reach 

is 

18000 

m. s^ = 

0.0005. 

Manning's roughness 

It 

coefficient = 0 

1 

o 

« 

Breadth 

of channel 

= 50 m 

1 . Assume the channel is 

rectangular in 

cross-section. 





Set # 4 

(Raudikivi, 

1979) 






The 

hydrograph 

at the 

upstream end of a 

r iver 

reach (x=0) 

J 


time t=0 

is as 

fol lows: 






t hours: 

0 

1 

2 

3 4 

5 

6 

7 8 

/ 

9 

Q m^/s : 

10 

20 

52 

100 141 

150 

144 

130 114 

100 
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t hours: 

10 

11 

12 

13 

14 

15 16 

17 

Q m^/s : 

86 

73 

61 

50 

41 

32 

25 19 

t hours: 

18 

19 

20 

21 

22 

23 24 

25 

Q m^/s : 

14 

11 

10 

10 

10 

10 

10 10 

t hours: 

26 

27 

23 





Q m^/s : 

10 

10 

10. 





Length of 

reach 

is 

18000 m. 

s ' = 

0.0005. 

Manning'^s 

roughness 


coefficient = 0.02. Breadth of Channel is 50 m. Assume a 

rectangular cross-section. 

Set # 5 

A rectangular channel 60.0 wide and 30 Km long carrying a 

^ 3 

steady uniform discharge of 50.0 m /s is subjected to upstream 
flood wave with a peak of 2000 m^/s increasing linearly in a period 

i 

of 10 hours. This upstream flow decreases linearly from its peak 

3 3 

of 2000 m /s/^50 m /s in 20 hours. Additional properties are s^ = 
0.0015 and n = 0.02. 

Set # 6 

The hydrograph at the upstream end of a river reach (x=0) 


is as follows: 
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T hours 

Dischargetm^/s) 


T hours 

Discharge Cm^/sl 

0 

10.00 


15 

60.24 

1 

25.28 


16 

43.48 

2 

65.28 


17 

28.28 

3 

114.72 


IS 

17.92 

4 

154.72 


19 

11.75 

5 

170.00 


20 

10.00 

6 

168.25 


21 

10.00 

7 

163.08 


22 

10.00 

8 

154.72 


23 

10.00 

9 

143.53 


24 

10.00 

10 

130.00 


25 

10.00 

11 

114.72 


26 

10.00 

12 

140.63 


27 

10.00 

13 

123.30 


28 

10.00 

14 

85.28 




Length of reach is 30000 

nri. 

s^ = 0.0005. Manning's roughness 

coefficient 

= 0.02. Breadth 

of Channel 

is 50 m. Assume a 

rectangular 

cross-section. 



• 

Set # 7 





The hydrograph at the 

upstream end of 

a river reach (><=0) 

is as follows; 


- 


T hours 

Discharge (m^/s) 


T hours 

Discharge (m^/s) 

0 

10.00 


15 

50.00 

1 

25. 28 

i 

16 

36.47 



100 


2 

65.28 

17 

53.78 

3 

114.72 

18 

70.90 

'4 

154.72 

19 

95.43 

5 

170.00 

20 

65. 75 

6 

168.25 

21 

51.00 

7 

163.08 

22 

30.60 

8 

154.72 

23 

15.23 

9 

143.53 

24 

10.00 

10 

130.00 

25 

10.00 

11 

114.72 

26 

10.00 

12 

98.36 

27 

10.00 

13 

81.64 

28 

10.00 

14 

65.28 



Length of reach is 30000 m. 

s^ = 0.0005. Mann Inga's roughness 

coef f icient 

= 0.02. Breadth 

of Channel 

is 50 m. Assume a 

rectangular 

cross-section. 



Set # 3 




The hydrograph at the upstream end of 

a river reach tx=0) 

is as follows: 



T hours 

Discharge (m /s) 

T hours 

Discharge (m^/s) 

0 

10.00 

15 

50.00 

1 

25.28 

16 

36.47 

2 

65. 28 

17 

- 43.78 

3 

114.72 

18 

65.90 

4 

154.72 

19 

81.43 



5 

170.00 

h 

168.25 

7 

163.08 

8 

154.72 

9 

143.53 

10 

130.00 

11 

114.72 

12 

98.36 

13 

81.64 

14 

65.28 


Length of reach is' 30000 m, s^ 
coefficient = 0-02. Breadth of 
rectangular cross-section. 
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20 

55.75 

21 

31.00 

22 

20.60 

23 

10.00 

24 

10.00 

25 

10.00 

26 

10.00 

27 

10.00 

28 

10. 00 


= 0.0005'. Manning''s roughness 

Channel is 50 m. Assume a 



